{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e964989e",
   "metadata": {},
   "source": [
    "# Integration\n",
    "**강좌**: *수치해석*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98b0dd7b",
   "metadata": {},
   "source": [
    "## Newton-Cotes\n",
    "### 사다리꼴 공식\n",
    "- 적분 구간을 사다리꼴 면적으로 계산\n",
    "    \n",
    "    $$\n",
    "    I = \\int_a^b f(x) dx \\approx \\int_a^b f_1(x)dx \n",
    "    = \\int_a^b \\left [\n",
    "    f(a) + \\frac{f(b) - f(a)}{b -a} (x-a) \n",
    "    \\right ] dx\n",
    "    $$\n",
    "    \n",
    ":::{figure-md} Trapezoid\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/4/40/Trapezoidal_rule_illustration.svg\" alt=\"trapezoid-fig\">\n",
    "\n",
    "Trapezoid rule (from Wikipedia)\n",
    ":::    \n",
    "\n",
    "- 오차\n",
    "    \n",
    "    $$\n",
    "    E_t = -\\frac{1}{12} f''(\\xi) (b-a)^3\n",
    "    $$\n",
    "    \n",
    "    * 유도\n",
    "        - Newton 보간식\n",
    "            - $\\alpha = \\frac{x-a}{h}$, $h=b-a$, $\\Delta f(a) = f(b) - f(a)$ 일때\n",
    "            \n",
    "            $$\\begin{align}\n",
    "            I &= \\int_a^b \\left [\n",
    "            f(a) + \\Delta f(a)\\alpha + \\frac{f''(\\xi)}{2} \\alpha (\\alpha-1)h^2\n",
    "            \\right ] dx\n",
    "            \\end{align}\n",
    "            $$\n",
    "        \n",
    "            - $f''(\\xi)$ 가 거의 변하지 않는 상수로 가정하면\n",
    "        \n",
    "            $$\\begin{align}                \n",
    "            I & \\approx h \\left [\n",
    "            \\alpha f(a) + \\frac{\\alpha^2}{2} \\Delta f(a) + \\left (\n",
    "            \\frac{\\alpha^3}{6} - \\frac{\\alpha^2}{4}\n",
    "            \\right) f''(\\xi) h^2\n",
    "            \\right ]_0^1 \\\\\n",
    "            & = h \\left [\n",
    "            f(a) + \\frac{\\Delta f(a) }{2} \n",
    "            \\right ] - \\frac{1}{12} f''(\\xi) h^3 \\\\\n",
    "            & = h \\frac{f(a) + f(b)}{2} - \\frac{1}{12} f''(\\xi) h^3 .\n",
    "            \\end{align}\n",
    "            $$\n",
    "        \n",
    "### 예제\n",
    "$f(x) = 0.2 + 25x - 200 x^2 + 675x^3 -900 x^4 + 400 x^5$ 을 $[0, 0.8]$ 구간에서 적분하시오.\n",
    "\n",
    "- Exact Integration $I_{exact} = 1.640533$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "90a0e2f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.dpi'] = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1cc63796",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1728000000000225\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAJqCAYAAAASQ256AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAACAjklEQVR4nO3deXxU1f3/8feZrGQjLCEgsm8BQVAENxREa93qWnGru22t2oq2frvvtfpt+9Pa9lutRSuKuOGuSN1QwQUQAZGdACEQlgQSEgIJSeb8/hgymSEBEnKTM8vr+Xj44J4zM3fecwzhfuaee66x1loBAAAAQJTyuQ4AAAAAAK1BUQMAAAAgqlHUAAAAAIhqFDUAAAAAohpFDQAAAICoRlEDAAAAIKpR1AAAAACIahQ1AAAAAKIaRQ0AAACAqEZRAwAAACCqUdQAAAAAiGoUNQAAAACiGkUNAAAAgKhGUQMAAAAgqiW6DuDa1q1bZa119v5du3aVJJWUlDjLEEsYT+8xpt5iPL3HmHqL8fQeY+otxtN7kTCmxhh17979iF8f90WNtdZpUROaA95hPL3HmHqL8fQeY+otxtN7jKm3GE/vRfOYMv0MAAAAQFSjqAEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENUoagAAAABENYoaAAAAAFGNogYAAABAVKOoAQAAABDVKGoAAAAARDWKGgAAAABRjaIGAAAAQFSjqAEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENUoagAAAABENYoaAAAAAFEt0XUAAAAAr9ldpVLhetk9u2V8Pqn+P5PQsO3zSUlJ0lF9ZNLSXUcG0AqeFTVvvPGGVq5cqY0bN2rXrl2qqalRdna2hg0bposuuki9evVq9r5uv/12FRcXH/TxBx98UD179vQiNgAAiGLWWql4q1S4TnZj4D8VrpN2lTY853A7MUY6qrdM/yHSgDyZAXlSbk8ZY9o0OwDveFbUvPzyy6qqqlKfPn3Uu3dvSVJhYaE++ugjffLJJ7rnnnt03HHHtWif48ePb7I/LS2t1XkBAEB0suVlsp99ILtkfqCA2bunlTu00uYC2c0F0py3A0VQeqbUf4hM/yEyw0ZJ/QZT5AARzLOi5p577lH//v2VnJwc1v/2229rypQpeuSRR/Twww/L52v+ZTy33367V/EAAEAUs3V10rIv5P/4XWnJfKmu7vAv6thZ6tQlULT46yS/P/Cf9Tds79kt7als/NrKCmnp57JLP5d99Wmp7yCZr10kc/wpMonM3gcijWd/K/Py8prsP/vss/XGG29o69atKioq0tFHH+3VWwIAgBhntxfJzn1X9tP3pbKdB39itx4yvfpLvfvL9N7/Z1anw+/fWmnbZtn8lVL+Stl1q6SijYFCKNSGNbL//otsp64yZ14gc9rZMmkZrfx0ALzSLl811J+dSeSbDQAAcBjWXye7YK7sR/+VVn/V9JOyO8uccqbMMcdLvfrJdDiyqenGGKn70TLdj5ZOPSvw/nsqpfWrZfNXyq5YIq1d3vCC0hLZGU/Ivv6szKlnyZz5DZluPY7ovQF4p82rjA8//FBFRUXq0aOHunXr1qLXvvbaa9q6dauSkpLUq1cvjR07VllZWW2UFAAAuGZXLZX/uSlS4frGDyYkSCPHyjfua9Ixx8n4Etokg0lLD+z/mOOkC6+S3Zgv+85rsgvmSHW1gSdVV8m+/4bs7DelkSfKd/E1Mj37tEkeAIdnrD3w/GrrvPbaayosLFR1dbU2b96swsJCderUST/+8Y/Vv3//Zu3jYKufpaSk6MYbb9TEiRObnefuu+9u1JecnKz7779fklRTU9PsfbWF+rNXtbW1TnPECsbTe4yptxhP7zGm3nI1nrVbNqli6j9U/dmHjTP16qcOZ31DHSacI1/Hw08payt1O4u1Z+aL2vPfl2UrysMfTExUxhU3K/3Sb8kkJB7wED+jXmI8vRcpY5qUlHTEr/X8TM2SJUu0dOnSYLtLly76/ve/3+yCRpJGjx6t4cOHq3///srKytK2bds0e/ZszZw5U4888ogyMjI0duxYr6MDAIB25q/crd0v/Ed73nhBqg3/ojH1tK8p7RtXKGnQsIhYeSyhc44yv3WrMi6/QXtnz1Tla8+prmhj4MHaWu1++l+qmveRsu/8pRJ79XMbFogznp+pqVdZWamNGzdqxowZWrp0qa688kpdeumlrdrnu+++q0cffVQ9evTQQw895EnOLVu2qI2GoFlycnIk6ZD35UHzMZ7eY0y9xXh6jzH1VnuNp62rk53ztuxr06WKXeEP9h8i36SbA/eLiWDW75cWz5P/2X9LpSUNDyQmylx4jczZF8skJPAz6jHG03uRMKbGGPXoceTXpzV/feUWSk9P19ChQ/XTn/5U/fv313PPPae1a9e2ap8TJ05Ux44dtWXLFm3fvt2jpAAAoD3ZNcvl//1k2acfDi9oOneVueWH8v3kTxFf0EiS8flkjj9Zvt/8XWb/IgOSpNpa2Zemyv+/P5bdUugsHxBP2qyoqZeYmKhTTjlF1lotXLiwVfvy+XzKzc2VJJWWlh7m2QAAIJJYf538bzwn/59/Jm0uaHggJVXmomvk+93D8p04PiKmmrWESUuX74YfyPeDX0nZnRseWL9a/t9NVuXL0wL32QHQZtpljeXMzExJUnl5+WGeeXiVlYEbZKWmprZ6XwAAoH3Ysp3yP/aAtPLLsH5zypkyl3xLJruLo2TeMSNOkO83/5B97t+yn84OdNbWqGLq/6lq3key375HJivbZUQgZrX5mRpJWr48sL57/VmWI1VYWKiioiKlpKSoZ8+eXkQDAABtzC5bJP/v7gwvaLrmyveTP8l3450xUdDUM+kZ8t10l3x3/EIKWamtZuVS+e//H9ntWxymA2KXJ0XNihUr9Mknn6jugFOrtbW1euutt/TRRx8pOTlZp5xySvCxWbNmafLkyZo+fXrYaxYvXqx169Y1eo+CggI98MADstZq4sSJ3MgTAIAIZ2tr5X9pqvx//XXYtTNm9Kny/fKvUXHdzJEyI8fK99t/yJw4vqGzeGugsCnIdxcMiFGeVAbbtm3TP//5T2VmZqp///7KzMxURUWFNm7cqNLSUiUlJem2225T165dg68pLy9XUVFRo2tjVq9erRkzZignJ0e5ubnKysrS9u3btX79etXV1WnYsGG6+uqrvYgNAADaiN1RLP+//yzlr2zoTEySufLbMqd/PequmzkSJj1TuvlupQ85Rruf/Gegs2KX/H/+mXy3/VRm2Cin+YBY4klRM2zYMF1yySVavny5Nm7cqPLyciUmJqpbt2468cQTdd5556l79+7N2teoUaO0Y8cO5efnq6CgQHv27FGHDh2Ul5encePG6YwzzpDP1y6z5gAAwBGwiz+T/z9/k/bsbujsfrR8371H5uj4un+LMUYZl16rhOzO2vWPP0p+v1S9V/6//U7mpsnyjT3ddUQgJnhS1HTr1k1XXXVVi14zadIkTZo0qVH/4MGDNXjwYC9iAQCAdmStlf3vS7IvTg3rNydPlLn6uzKpHRwlc6/DxPNVoQT5H7lf2lct1dXK/vsv8peXyXfWha7jAVGPUx4AAKDVrN8v+/xj4QVNSqrMTXfJd9PkuC5o6pkRo+X74R+kjMxgn31uivwvTnV6I3AgFlDUAACAVrE1NbJT/p/su681dOZ0l+/nD8h38hnugkUg03+IfD/+X6lLt2CfnfWi7H/+Kltb6zAZEN0oagAAwBGze/fI/7ffyi6Y09DZe4B8P/lfmR5HuwsWwUz3o+X7yf9KPfsE++yns2Ufe0DW73eYDIheFDUAAOCI2F2l8v/5p+H3nxk2Sr577pXJ6nTwF0Imu4t8/3OfNHh4sM9+Plf2xSfchQKiGEUNAABoMbt1s/z3/49UuD7YZ04cL9/3fymTmuYwWfQwaRnyTf6NNOy4YJ99+xX533vDXSggSlHUAACAFrHr18j/vz+WSrYF+8zZF8vcdJdMYpLDZNHHJCXL970fS70alrq2z/1b9otPHaYCog9FDQAAaDa7fJH8f/mZtLs82Gcuv1G+y2+S4T5yR8Skpsn3g19JnXMCHdbKP+X/yYbeuBTAIfHbBwAANItdtVT+f9wbuM+KJCUkytzyQ/nOvsRtsBhgsrvI94NfSx3SAx01++T/x+9ltxW5DQZECYoaAABwWDZ/pfx//71Usy/QkZIq3w9+Kd+J490GiyGmZ2/5bv+ZlLj/3ui7K+R/6Dey5WVOcwHRgKIGAAAckt24Tv6//Vaqrgp0JCfL94NfyYRc4A5vmCEjZG64s6GjeKv8//iDbHW1u1BAFKCoAQAAB2WLNsr/4K+kPZWBjsRE+W77uUzIUsTwlu/E8TKXXt/QsX61/FP+IuuvcxcKiHAUNQAAoEl2+xb5H/hVw6IAPp983/0fmWM4Q9PWzDmXykw4t6Fj8TzZZ//tLA8Q6ShqAABAI3ZHsfwP/FLatTPQYUxgyeZRJ7kNFieMMTJXfkcaOTbYZ2fPlP/T2Q5TAZGLogYAAISpK90RKGh2bA/2mevuYFGAdmYSEuT79o+kvoOCffbpR2S3syIacCCKGgAAEOQv36XSX/9ACjlwNld+R75xX3OYKn6ZlFT5bv2JlLZ/qefqvfI/+hfZ2hq3wYAIQ1EDAAAkSbZqr3b+drJqN64L9plLr5PvzAscpoLpkiPfdd9v6ChYK/vKNHeBgAhEUQMAAGT9dfL/+y+qDbmLvTl/knznftNhKtQzo0+ROf2cYNv+92XZZYscJgIiC0UNAACQfeEJ6csFwbY58xsyF13jLhAaMZNulnr0Crb9jz8oW17qMBEQOShqAACIc/4PZ8m++2qwnXLi6TKTbpYxxmEqHMikpMj3nXukxKRAR3mZ/P95SNbvdxsMiAAUNQAAxDG7fJHs9EeC7cT+g9Xxrt/I+DhEiETm6L4yk25q6PjqC9l3X3MXCIgQ/MYCACBO2S2F8j/yJ6n+m/7szur08z/Ll9rBbTAckplwXvj9a156UrZgrbtAQASgqAEAIA7Zil3y/+130t7KQEdyinx3/FIJXbq5DYbDMsbId8MPpOwugY662sAyz1V73QYDHKKoAQAgztiaGvn/+UepZFugwxj5bvmhTJ8BboOh2UxGlny33C3VX/e0vUj2mUfdhgIcoqgBACCOWGtln/y7tHZFsM9cdr3McSc5TIUjYYaMkDn38mDbfvKe/PM/cpgIcIeiBgCAOGLffE72sw+CbXPqWTJnX+IuEFrFfONKaUBesG2feVS2ssJhIsANihoAAOKE/Xyu7KvTGzqGjJD51vdYujmKmcRE+W6+W0pOCXTsLpd9carbUIADFDUAAMQBW7RR/if+1tCR21O+7/1Epv6eJ4haJqe7zIVXBdt2ztuya5c7TAS0P4oaAABinK3aK/8j/ytVVwU6OqTL9/1fyqRnug0Gz5gzL5R69gm2/dMelq2tdZgIaF8UNQAAxLDAwgD/kLYUBvt8N98lk3uUw1TwmklMlO9btzV0bC6QfY+bciJ+UNQAABDD7AczZRfMCbbNuZfJhNy4EbHDDBwqc/rXg2372jOyO7Y7TAS0H4oaAABilF23Sva5xxo6hoyQuehb7gKhzZlLr5MyOwYa+6rln/4vWWvdhgLaAUUNAAAxyFaUy/+v/5Xq9l9X0bGzfN/+kUxCgttgaFMmPVPm8psaOr5cIC2e5y4Q0E4oagAAiDHW75f/sf8n7SwJdPh88n3nHpmOndwGQ7swJ02QhowItv3PPCpbtcddIKAdUNQAABBj7JvPS8sWBdvm0utlBh/jMBHakzFGvm99T0pMDHSUlsi+9ozbUEAbo6gBACCG2GWLZF8POYA97iSZsy92lgdumO5Hy5xzWbBt33tdduM6h4mAtkVRAwBAjLA7i+Wf8hep/sLwbj3ku+FOGWPcBoMT5rzLpZzugYbfL/+0f8r669yGAtoIRQ0AADHA1tbK/68/SbsrAh1JyfLd+hOZtHS3weCMSUqW75rvNXSsXy370dvuAgFtiKIGAIAYYF9/Vlq3Ktg213xPplc/h4kQCcwxx8mMOS3Yti89KVtR7jAR0DYoagAAiHJ29Veyb70QbJtTz5Lv1DMdJkIkMVfcInXYf8Zub6XszBcO/QIgClHUAAAQxeye3fI/9mDIdTRHyVz5bbehEFFMx04y508Ktu0Hb8qWbHOYCPAeRQ0AAFHKWis77WFpZ3GgIyFBvlt+KJPawW0wRBwz8XypU9dAo7ZW9tXpbgMBHqOoAQAgStl5H8gumBNsmwuvluk3yGEiRCqTlCxz0dXBtp33geym9Q4TAd6iqAEAIArZ4q2yTz/S0DH4GJlzLnUXCBHPnHyGdFTvQMNa+V96ym0gwEMUNQAARBlbVyf/Yw9IVXsDHR3S5bvpbhlfgttgiGjGlyDfpdc1dCz9XHbVUneBAA9R1AAAEGXszBek/JXBtrn2NpkuOQ4TIWocO0YaOCzY9L84VbZ+kQkgilHUAAAQRWz+Stk3ng22zclnyBdyHxLgUIwx8l12fUPH+tXSF5+6CwR4hKIGAIAoYffukX/K/5P8/kBH11yZq77rNhSijhk4VBp1YrDtf/kp2bo6h4mA1qOoAQAgSthnHpXq7y/i8wWWb+6Q5jYUopLvkmsls/8wcNtm2Y/fcRsIaCWKGgAAooB/wVzZT98Pts35V8gMyHOYCNHMHNVb5tQzg2372rOy1VUOEwGtQ1EDAECEs+WlstMfbugYkBd2h3jgSJhvXCUlJQcau3bKvvua20BAK1DUAAAQway18j/9iLS7ItCRkirfzXfLJLB8M1rHdO4qc+Y3gm3735dkK8odJgKOHEUNAAARzC6YE7Y6lbnsBpmc7g4TIZaYcy6T0tIDjb17AsuFA1GIogYAgAhly0tln/lXQ8eQETLjz3EXCDHHpGfInHd5sG0/eFN2x3aHiYAjQ1EDAEAEstbKP+3h8Gln139fxsc/3fCWOeN8qVPXQKO2Vvb1Zw/9AiAC8ZsRAIAIZBfMkRZ9Fmwz7QxtxSSnyFx4VbBtP/tAdmexw0RAy1HUAAAQYZh2hvZmTjpD6pwTaNTVyr79itM8QEtR1AAAEEGYdgYXTGKizDmXBtt2zn9ly8vcBQJaiN+QAABEEKadwRVz6llSZsdAY98+2fdedxsIaAGKGgAAIgTTzuCSSU6R+drFwbad/absnkp3gYAWSPRqR2+88YZWrlypjRs3ateuXaqpqVF2draGDRumiy66SL169WrR/iorK/XCCy9o/vz5KisrU3Z2tsaMGaNJkyYpPT3dq9gAAEQEpp0hEpgJ58q+NUPaWxm4b80HM8OWfAYilWe/KV9++WUtWrRIGRkZGjFihI4//nglJSXpo48+0o9//GMtWrSo2fuqqKjQz372M82cOVMJCQkaM2aMOnTooLfeeks//elPVVFR4VVsAAAigp3/EdPO4JzpkCYz8fxg2777mmx1tcNEQPN4dqbmnnvuUf/+/ZWcnBzW//bbb2vKlCl65JFH9PDDD8vXjG+cpk6dqi1btmjs2LG66667lJCQIEl6/PHHNWvWLE2dOlV33HGHV9EBAHDKlpfJPvtoQwfTzuCQOfNC2XdelfZVSxW7ZOe+I3PmBa5jAYfk2ZmavLy8RgWNJJ199tnq3r27SktLVVRUdNj9lJWVac6cOUpISNAtt9wSLGgk6dprr1VWVpbmzp2rsrIyr6IDAOCUff4xpp0hYpjMLJnTvx5s2/++JFtb4zARcHjt8huz/uxMYuLhTwwtWrRI1loNGzZM2dnZYY8lJSVp9OjR8vv9Wrx4cRskBQCgfdmvvpCd92GwbS65lmlncM6cfYmUsP+4rbRE9rMPnOYBDqfNi5oPP/xQRUVF6tGjh7p163bY5xcUFEiS+vXr1+Tj9f0bNmzwLCMAAC7Y6ir5p/2zoaPfYJkzznMXCNjPdOoic8rEYNu+9aKsv85hIuDQPLumpt5rr72mwsJCVVdXa/PmzSosLFSnTp105513Nut6mpKSEklS586dm3y8S5cuYc87nLvvvrtRX3Jysu6//35JUteuXZu1n7ZSf/YqJyfHaY5YwXh6jzH1FuPpvWge0/In/q49O7YHGgkJ6nLnL5WU6/YsTTSPZ6SK1jGtvfoWlXz8ruT3S9uLlLnmK3UYd5brWFE7npEsFsbU86JmyZIlWrp0abDdpUsXff/731f//v2b9fqqqipJUkpKSpOP1/dXsxIHACCK1eSv0p7Xng220y++Rkl9BzpMBIRL7NFLqaeeqao570iSKmc8qdRTz5QxxnEyoDHPi5pf/vKXkgL3mdm4caNmzJih3/zmN7ryyit16aWXHvb11lpJ8uwvzAMPPHDIx0tKSoLv6UJ9RVxcXOwsQyxhPL3HmHqL8fReNI6prauT/6HfB74Bl6Sc7to78RuqioDPEI3jGemieUztxAuk/UVN7YY1Kp49S2bECU4zRfN4RqpIGFNjjHr06HHEr2+za2rS09M1dOhQ/fSnP1X//v313HPPae3atYd9XYcOHSQ1nLE5UP0ZmoOdyQEAINLZ916XNuYH275rb5dJ5t81RB5zdD9p5Nhg2//m806/DAYOps0XCkhMTNQpp5wia60WLlx42OfXX+Oyc+fOJh/fsWNH2PMAAIgmtnir7KtPB9vmlDNlho50mAg4NN+532xo5K+UVi9zFwY4iHZZ0jkzM1OSVF5eftjn9unTR5K0fv36Jh+v769/HgAA0cJaK//0RwI3NZSkjCyZy290Gwo4DDMgTxoyItj2z3zeYRqgae1S1CxfvlySlJube9jnjho1SsYYrVixQrt27Qp7rKamRgsXLpQxRscdd1ybZAUAoK3Y+R9JX30RbJsrbpHJyHKYCGge33mXNzSWL5bd1PSXz4ArnhQ1K1as0CeffKK6uvD1y2tra/XWW2/po48+UnJysk455ZTgY7NmzdLkyZM1ffr0sNd06tRJp556qmprazVlypSwfU6bNk3l5eUaN25coxtzAgAQyezuctnnpjR0HHOczInj3QUCWmLoSKn3gGDTvv+mwzBAY56sfrZt2zb985//VGZmpvr376/MzExVVFRo48aNKi0tVVJSkm677baw62DKy8tVVFSk0tLSRvu74YYbtGbNGs2bN0+TJ0/WgAEDVFhYqMLCQuXm5ur666/3IjYAAO3GzviPVLF/BkJysnzXfI+lcRE1jDEyZ35D9j9/lSTZeR/IXnodZxoRMTwpaoYNG6ZLLrlEy5cv18aNG1VeXq7ExER169ZNJ554os477zx17978m4llZWXpvvvu0/PPP68FCxZo/vz56tixo8455xxNmjRJGRkZXsQGAKBd2JVfyn78XrBtLrxGJsftTTaBljJjTmsozvftk537jsw5l7mOBUjyqKjp1q2brrrqqha9ZtKkSZo0adJBH8/IyNBNN92km266qbXxAABwxtbWyP/0ww0dvfvLnHWhu0DAETJJSTKnf132zcBCAXb2TNmvXSyTkOA4GdBOCwUAABCv7NuvSFs3BxrGBO5Jw0EgopQZf65U//O7s1haMs9tIGA/ihoAANqILdkm++ZzwbaZcK5M30HuAgGtZDp1kTm+YeEn/3tvOEwDNKCoAQCgjfif/be0b1+gkdlR5uJvuQ0EeMCc+Y2GxuqvWN4ZEYGiBgCANmAXz5OWzA+2zaSbZNJY6AYxoP8Qqc/AYJPlnREJKGoAAPCYra4OnKWpN3i4zIkTnOUBvFS/vHM9+9kHsrvLHSYCKGoAAPCcnfm8tGN7oJGQIN81t3JPGsQUc8I4KbNjoFGzT3bOO24DIe5R1AAA4CG7pVD2vy8H2+ZrF8sc1dthIsB7JilJZvw5wbb94E3ZujqHiRDvKGoAAPCItVb+px+R6moDHZ27ylxwhdtQQBsx488JWd65RFrM8s5wh6IGAACP2PkfSauWBtu+K78jk5LqMBHQdkx2F5nRpwbb/vdfd5gG8Y6iBgAAD9g9lbLPP9bQMeIEadSJ7gIB7cBMvKChsXqZ7MZ17sIgrlHUAADgAfvq01J5WaCRlCzfVd9hcQDEvv5DpJAbytr3uRkn3KCoAQCglWxBvuzsmcG2OX+STE53h4mA9hFY3rnhbI2d/5FsBcs7o/1R1AAA0ArW75f/6Ycl6w905PaUOfsSt6GAdmRGj5OysgONmn2yc992mgfxiaIGAIBWsJ+8J61fHWz7rv6uTFKSw0RA+2q0vPPsmSzvjHZHUQMAwBGye3bLvvRksG1GnyozbJS7QIAj5vSQ5Z1LS6TFn7kNhLhDUQMAwBGyrz0jVewKNJJTZCbd5DYQ4IjJ7hyYhraf/6P/OkyDeERRAwDAEbCbNsjOfjPYNuddLtM5x2EiwK3QKWhavli2eKu7MIg7FDUAALSQtVb+Zx6V/PsXB+jWg8UBgEHDpO49g037yXsOwyDeUNQAANBCdsEcafVXwbbviltYHABxzxgjM+7sYNvOfVfWz4IBaB8UNQAAtICt2iv7wn8aOo4dI3PsGHeBgAhiTj6jYcGAsh3SskVuAyFuUNQAANACdubzgYM1SUpMlO+Km90GAiKIycqWRp4YbPvncM8atA+KGgAAmslu3Sz79qvBtjn7UpluRzlMBEQe32lfa2h8uUB2V6m7MIgbFDUAADSDtVb+56ZIdbWBjs5dZc77pttQQCQaNkrq3DWwXVcn++n7TuMgPlDUAADQHF8ukL5aGGz6Lr9JJiXVYSAgMhlfgsypZwXbds47stY6TIR4QFEDAMBh2Jp98j/774aOISOk0ae6CwREOHPqWZIxgcb2Imn1MreBEPMoagAAOAz735ekkm2Bhs8n31Xflak/YAPQiOnSLTANbT87lwUD0LYoagAAOAS7Y7vsWzOCbTPxGzI9eztMBEQH32kh96xZ+Ils5W6HaRDrKGoAADgE+8J/pH37Ao2sbJlvXOk2EBAtRo6VMrIC2zX7ZOd/6DYPYhpFDQAAB2FXfSW78ONg21x6vUxausNEQPQwiUkyp0wMtu1Hb7NgANoMRQ0AAE2w/rrwxQH6DgrcLR1As5lxIfes2bRe2pjvLgxiGkUNAABNsHPeCRyE7ee78tsyPv7ZBFrC9OglDRwabNs5LBiAtsFvZwAADmArd8u+8lSwbU46Q2ZAnsNEQPQy40IWDJj/kWx1lcM0iFUUNQAAHMC+/oy0uyLQSEmVuew6t4GAKGZOOFVK7RBo7N0Tdp0a4BWKGgAAQtiijbKz3wy2zXmXy2R3cZgIiG4mJVVm7Phg2855x2EaxCqKGgAA9rPWyv/cFMnvD3TkdJf52kVuQwExwJwWsmDA2uWyWza5C4OYRFEDAEC9JfOl5YuDTd/lN8kkJbvLA8SKPgOlo/sFm3YuZ2vgLYoaAAAk2Zoa+Z9/rKFj6Ehp1InuAgExxBgTdrbGfvq+bG2tw0SINRQ1AABIsu+9JhVvDTR8PvmuuEXGGLehgBhiTpwgJSYFGhW7pBWLXcZBjKGoAQDEPVu2U/aN54NtM/5cmZ59HCYCYo9Jz5COHRNs288+cBcGMYeiBgAQ9+zLT0nVewON9EyZi652GwiIUb6TJgS37eLPZKv2uAuDmEJRAwCIa3b9atlP3gu2zUXXyKRnOkwExLDho6W0jMD2vn2yX3zmNg9iBkUNACBuWWvlf/bfDR09+8ic/nV3gYAYZ5KSZE4YF2zbeR+4C4OYQlEDAIhbdt6H0rpVwbbvym/LJCQ4TATEPhMyBU0rvpQt2+EsC2IHRQ0AIC7Z6irZF6c2dBx/skzese4CAfFiQJ7UpVtg2/pl53/kNg9iAkUNACAu2VkvSfXfECcmyvfNG90GAuKE8fkCyzvvxypo8AJFDQAg7tgdxbL/fSnYNmddJJPT3WEiIL6Yk8Y3NArXy24ucBcGMYGiBgAQd+xLU6WafYFGVrbM+Ze7DQTEGdOjl9RnYLDNggFoLYoaAEBcsWtXhM3hN5dcK5Oa5jAREJ9Cz9bYeR/K+v0O0yDaUdQAAOKG9fvDl3DuPUDmlDPdBQLimBlzumT2H4ruLJHWLHcbCFGNogYAEDfsZ7OlgrXBtu+KW2R8/FMIuGA6dpKOGRVsMwUNrcFvcgBAXLBVe2VfeirYNqNPlRl8jMNEAMJWQfv8Y9n6a92AFqKoAQDEBfvWi9KunYFGYpLMN29wmgeAZI47SUpJDTT2Vkpffu42EKIWRQ0AIObZkm2yb78cbJuzL5HpmuswEQBJMimpgcJmPz/3rMERoqgBAMQ8O+MJqbYm0OjYWebcy5zmAdAgdAqaln4uW1nhLAuiF0UNACCm2dXLZBd+HGybS6+VSe3gMBGAMENHSlnZge26WtnPPz7k04GmUNQAAGKW9fvlf25KQ0ffQTInneEuEIBGTEKCzNjTg23LFDQcAYoaAEDMsp+8J23MD7ZZwhmITOakCQ2Ntctli7c6y4LoxG92AEBMslV7ZF8OWcJ5zGkyA4c6TATgoHoPkLofHWza+R85DINolNjaHVRXV2vJkiVauHCh8vPzVVxcLL/fr+7du+vEE0/UBRdcoNTU1Gbv7/bbb1dxcfFBH3/wwQfVs2fP1sYGAMQ4O3OGVF4WaCQny1x2g8s4AA7BGCNz0gTZV6ZJCkxBs+ddLmOM42SIFq0uaubOnat//etfkqRevXpp5MiR2rt3r1avXq3nn39eH3/8sX7zm9+oY8eOLdrv+PHjm+xPS0trbWQAQIyzxVtl33k12DZnXyrTJcdhIgCHY8aeHixqtHWTVLBW6jvIbShEjVYXNYmJiTr77LN1/vnnq0ePHsH+0tJS3X///Vq/fr2eeOIJ3XnnnS3a7+23397aaACAOOV/8YmGJZyzO8ucc6nTPAAOz+R0lwYOk9YulyTZBXNkKGrQTK2+pmb8+PG65ZZbwgoaSerUqZNuvvlmSdL8+fNVW1vb2rcCAOCw7OqvpIWfBNvm0utlUpo/DRqAO2bsacFt+/nHstY6TINo0qYLBfTp00eSVFNTo4oKbqQEAGhb1l8XvoRzv8EyJzY9nRlA5DHHnyLVX0ezs1hav9ptIESNVk8/O5Rt27ZJkhISEpSRkdGi17722mvaunWrkpKS1KtXL40dO1ZZWVltERMAECPsJ+9LG9cF2yzhDEQX07GTNOgYafVXkiT7+VyZ/kMcp0I0aNOiZubMmZKkUaNGKSkpqUWvnTZtWlh76tSpuvHGGzVx4sQW7efuu+9u1JecnKz7779fktS1a9cW7c9riYmB/wU5OVzA6gXG03uMqbcYT+/Vj2mXjDSVvPq06ierpJ72NWWfdNrBX4gm8TPqPca0ZfZMOEfl+4sas+gzdb3tx2GroDGe3ouFMW2zouaLL77Q7NmzlZCQoCuuuKLZrxs9erSGDx+u/v37KysrS9u2bdPs2bM1c+ZMPfLII8rIyNDYsWPbKjYAIEpVznhS/rKdgUZyijKvu81tIABHJOXkCdK//59krfwl21SzepmShwx3HQsRztg2uAJr06ZN+uUvf6nKykrdcMMNOu+881q9z3fffVePPvqoevTooYceesiDlAFbtmxxehFafUV8qHvzoPkYT+8xpt5iPL2Xk5Oj2m1FKrn9yuCKZ+aCK+W76GrHyaITP6PeY0xbru4vP5dWLZUkma9dJN+km4OPMZ7ei4QxNcY0WnisJTyfaLxjxw798Y9/VGVlpS644AJPChpJmjhxojp27KgtW7Zo+/btnuwTABAbKqb+gyWcgRhiRp8a3LYLP5b1+x2mQTTwtKgpLy/XH/7wB5WUlGjChAm69tprPdu3z+dTbm6upMA9cAAAkKR9yxap+pPZwTZLOAPRzxx/smT2H6buLGEVNByWZ0XN3r17dd9992nz5s0aO3asbr311rCLurxQWVkpSUpN5R8rAEBgCefyx/7a0MESzkBMMB07SYOPCbbt5x87TINo4ElRU1NToz/96U/Kz8/XyJEjNXnyZPk8XkKzsLBQRUVFSklJUc+ePT3dNwAgOtlP3lftuoZvcFnCGYgd5oSQKWhfMAUNh9bq3/x+v18PPfSQli1bpqFDh+pHP/pRcFm4g5k1a5YmT56s6dOnh/UvXrxY69ata/T8goICPfDAA7LWauLEiYfdPwAg9tmqPbIvPxVsm7GnywzIc5gIgJeYgoaWaHV1MGvWLM2fP1+SlJmZqSlTpjT5vGuvvTZ488zy8nIVFRU1ujZm9erVmjFjhnJycpSbm6usrCxt375d69evV11dnYYNG6arr2Y1GwCAZGe+IJWXBRrJyTKXXu80DwBvmaz9U9D2r4JmP/+YLy5wUK0uanbv3h3cri9umnL55ZcHi5qDGTVqlHbs2KH8/HwVFBRoz5496tChg/Ly8jRu3DidccYZnk9rAwBEH1u8VfadV4Pt9Eu+paou0XvTOABNMyeMk60vahZ+LHv5jY4TIVK1uqiZNGmSJk2a5MlrBg8erMGDB7c2EgAgxvln/EeqrZUk+brkBIqait2HeRWAaGOOP1l2+r8k65dK909B278aLhCK0x4AgKhiVy2Vvvg02M687nb5Ujs4TASgrZisbGnI8GDbfj7XXRhENIoaAEDUsP46+Z8NuXaz32ClnvY1d4EAtLnwG3F+wipoaBJFDQAgati570qb1gfbviu/zRLOQIwLWwWttEQ1q79yGwgRiX8JAABRwe6plH1lWrBtTpog03+Iw0QA2sOBU9CqPn7fXRhELIoaAEBUsG8+L1XsCjSSU1jCGYgj5oRxwe2qT95nChoaoagBAEQ8u61I9r3Xg21z7mUynbo4TASgPZnjTgpOQfPvKGYKGhqhqAEARDz/C49LdYElnNW5q8zXLnEbCEC7MlnZUt6IYJspaDgQRQ0AIKLZ5YulJQ03dzaX3SCTkuIuEAAnQldBYwoaDkRRAwCIWLauTv7nH2voGDhUZsxp7gIBcCZ0FTT/jmJp3SrHiRBJKGoAABHLzvmvtLkg2PZdcYuMMQ4TAXDFZHYMm4LGjTgRiqIGABCRbOVu2VefDrbNKWfK9B3kMBEA18wJB9yI01qHaRBJKGoAABHJvv6MtLsi0EhJlbnkWreBADhnjjtZqr/hbtkOqWCt20CIGBQ1AICIY4s2ys5+M9g2535TJruzw0QAIoHJ7KikvGODbbtonsM0iCQUNQCAiGKtlf+5KVL9ykY53WXOvthpJgCRI/XE04PbdvFnDpMgklDUAAAiy5L50vLFwabv8ptkkpLd5QEQUVJCihoVbZTdVuQuDCIGRQ0AIGLYmprwJZyHjpRGneguEICIk9i9pxL7Dgy2OVsDiaIGABBB7LuvScVbAw2fjyWcATQp9GyNXURRA4oaAECEsGU7Zd98Ptg248+V6dnHYSIAkSr0uhqtWyW7q9RdGEQEihoAQESwLz0pVe8NNNIzZS662m0gABErsd9gqXNOoGGt7JL5bgPBOYoaAIBzdt0q2U/fD7bNRdfIpGc6TAQgkhljZI47Kdi2i1naOd5R1AAAnLJ+v/zP/ruho2cfmdO/7i4QgKhgQhcRWbFYtmqPuzBwjqIGAOCU/ewDaf3qYNt35bdlEhLcBQIQHQYdI6VlBLZra6WvvnCbB05R1AAAnLFVewLX0tQ7/hSZkLuFA8DBmIQEmZFjgm27iClo8YyiBgDgjJ35grRrZ6CRmCTf5Te6DQQgqphRIdfVLP1ctrbGYRq4RFEDAHDCbt8i+86rwbb5+iUyXXMdJgIQdY45TkpKDmzvrZRWfeU2D5yhqAEAOOF//rHAPHhJyu4ic+433QYCEHVMSmqgsNnPLuZGnPGKogYA0O7s0oVSyH0lzDdvCBycAEALhU1BWzxP1u93mAauUNQAANqVrakJX8J54DCZsacf/AUAcAjm2DGS2X9IW7ZTKljrNhCcoKgBALQr++5r0vaiQMP45LvqOzLGuA0FIGqZzCxp0LBg2y5iClo8oqgBALQbW7pD9s3ngm0z/hyZ3v0dJgIQC8xxDTfitItZ2jkeUdQAANqNnfEfqboq0MjIlLn4GreBAMQEM7KhqNGWQtmtm9yFgRMUNQCAdmFXfyU7/6Ng21xyrUx6psNEAGKFyekuHd0v2OZsTfyhqAEAtDlbVyf/9H81dPQZKDPua+4CAYg5TEGLbxQ1AIA2Zz98S9pcEGz7rvqOjC/BYSIAsSZ0aWetWyVbttNdGLQ7ihoAQJuyFbtkX3062DannCkzIM9hIgAxqVc/qUu3wLa1sl/OP/TzEVMoagAAbcq+/JS0pzLQ6JAmc9l1bgMBiEnGGJnjQm7EydLOcYWiBgDQZuz6NbJz3wm2zYVXyWR1cpgIQCwLm4K24kvZvXvchUG7oqgBALQJ6/fL/8y/JGsDHT16yUw432kmADFu4FApY/+qinW1sl8tdJsH7YaiBgDQJuyn70vrVwfbvqu+I5OY6DARgFhnEhJkRoxp6Pjyc3dh0K4oagAAnrN7dsu+ODXYNqNPlRk60mEiAPHCHHtCcNt+tVDWX+cwDdoLRQ0AwHP2lWlSxa5AIzlF5vKb3AYCED+GHScl7F8yfne5tH6N2zxoFxQ1AABP2YK1sh+8FWyb8y6X6ZLjMBGAeGLS0qVBxwTb9ssFDtOgvVDUAAA8Y/118k97uGFxgO49Zb5+idtQAOKOGREyBY3rauICRQ0AwDN2zjvShoapHr6rb5VJTHKYCEA8MseGLBawab3szmJ3YdAuKGoAAJ6wFbtkX3oy2DZjTmNxAABOmO49pW5HBducrYl9FDUAAE/YF5+Q9uwONFI7yExicQAA7oStgraUoibWUdQAAFrNrl0u+/F7wba56BqZ7C4OEwGId2FT0FYska2udhcGbY6iBgDQKrauTv6nH2noOLqfzBnnuwsEAJI0aJiU2iGwXbNPWvWl2zxoUxQ1AIBWsbPfkDZtCLZ919wqU3+PCABwxCQmBe5Zsx9T0GIbRQ0A4IjZsh2yr04Pts2pZ8kMHOowEQA0CJ2CZr9cIFu/3DxiDkUNAOCI2ecfl6r2BhrpmTKX3eA0DwCEMiOOl4wJNHaWSJs3OM2DtkNRAwA4InbFEtkFc4Jtc+m1MplZDhMBQDiT1UnqOyjYZmnn2EVRAwBoMVtTI//0kMUB+g2WGXe2u0AAcBBhSzt/ucBhErQlihoAQIvZt1+Wtm4ONIxPvmu+J+PjnxQAkceMCFnaed0q2Ypyd2HQZvgXCADQInZbkewbzwXbZsK5Mn0GuAsEAIfSu7+U3Tmwba3ssoVu86BNUNQAAJrNWiv/tH9KtTWBjo6dZS7+lttQAHAIxhiZEQ1T0MR1NTGJogYA0Gz209nSyoYb2Pmu+rZMWrrDRABweGHX1Xz1hWxtrcM0aAsUNQCAZrEV5bIvPNbQcewY6fhT3AUCgObKGyklJgW291ZK+Svd5oHnKGoAAM1iX3hc2l0RaKSkynf1rTL1938AgAhmUjtIeSOCbVZBiz2Jrd1BdXW1lixZooULFyo/P1/FxcXy+/3q3r27TjzxRF1wwQVKTU1t0T4rKyv1wgsvaP78+SorK1N2drbGjBmjSZMmKT2daQ4A0N7siiWyn74fbJuLrpHpkuMwEQC0jBlxguxXX0iS7NLPpctvdJwIXmr1mZq5c+fqL3/5i2bPni1rrUaOHKm8vDxt375dzz//vH76059q165dzd5fRUWFfvazn2nmzJlKSEjQmDFj1KFDB7311lv66U9/qoqKitZGBgC0gK3ZJ/+0hxs6eg+QmXiBu0AAcATCFgvYUii7fYu7MPBcq8/UJCYm6uyzz9b555+vHj16BPtLS0t1//33a/369XriiSd05513Nmt/U6dO1ZYtWzR27FjdddddSkhIkCQ9/vjjmjVrlqZOnao77rijtbEBAM1k33xe2l4UaBiffNfdIbP/dzMARAuT0106qrdUtFFS4GyNOfMbjlPBK60+UzN+/HjdcsstYQWNJHXq1Ek333yzJGn+/PmqbcYqE2VlZZozZ44SEhJ0yy23BAsaSbr22muVlZWluXPnqqysrLWxAQDNYIs2ys56Kdg2Z36De9IAiFqhZ2ssSzvHlDZdKKBPnz6SpJqammZNG1u0aJGstRo2bJiys7PDHktKStLo0aPl9/u1ePHiNkgLAAhl/X75n/o/qW7/l1Kdc2QuutptKABoBXPsmIbG6qWyVXvchYGn2rSo2bZtmyQpISFBGRkZh31+QUGBJKlfv35NPl7fv2HDBm8CAgAOys59W1q7Itj2XXNrYAUhAIhWA/KktP3HpLW10vIlbvPAM62+puZQZs6cKUkaNWqUkpKSDvv8kpISSVLnzp2bfLxLly5hz2uOu+++u1FfcnKy7r//fklS165dm72vtpCYGPhfkJPDKkJeYDy9x5h6K1rGs25niUpefDLYTj1lorLPPM9hooOLljGNFoyn9xhTb7V2PMtGn6yqOe9IklJWL1XHr1/oWbZoFQs/o212puaLL77Q7NmzlZCQoCuuuKJZr6mqqpIkpaSkNPl4fX91dbU3IQEATap47K+ye3ZLkkxahjJvuctxIgDwRsoJpwa3q7/4TNZah2nglTY5U7Np0yb9/e9/l7VW1157rfr27dus19X/UHl5M7cHHnjgkI+XlJQ4/WGur4iLi4udZYgljKf3GFNvRcN42sWfyf/xew0dl1yrnXVWitDM0TCm0YTx9B5j6q3WjqftNVAyRrJW/tISFS9eIHN005c+xItI+Bk1xjRaeKwlPD9Ts2PHDv3xj39UZWWlLrjgAp13XvOnK3ToEJirXX/G5kD1Z2gOdiYHANA6tnK3/NMeaegYkCdz+tfdBQIAj5nMLKnvoGC7/oaciG6eFjXl5eX6wx/+oJKSEk2YMEHXXntti15ff33Lzp07m3x8x44dYc8DAHjLvvCYtGv/7+DEJPmu/4GMr03XlAGAdmeOOT64TVETGzz7l2rv3r267777tHnzZo0dO1a33npri6eR1S8BvX79+iYfr++vfx4AwDv2q4WyIdPOzIVXy/Q42mEiAGgbZnhDUaO1K1jaOQZ4UtTU1NToT3/6k/Lz8zVy5EhNnjxZviP4Zm/UqFEyxmjFihXatWtXo/dYuHChjDE67rjjvIgNANjP7t0TuCdNvT4DZc6+2FkeAGhTfQc1LO1cVyutXOo2D1qt1UWN3+/XQw89pGXLlmno0KH60Y9+FFwW7mBmzZqlyZMna/r06WH9nTp10qmnnqra2lpNmTJFdXV1wcemTZum8vJyjRs3rtGNOQEArWNffELauX+5/IRE+W74gUxCgtNMANBWTEKCzLBRwbZdxhS0aNfq1c9mzZql+fPnS5IyMzM1ZcqUJp937bXXKisrS1Lg2puioiKVlpY2et4NN9ygNWvWaN68eZo8ebIGDBigwsJCFRYWKjc3V9dff31rIwMAQtiVX8p+OCvYNudPkjm6r7tAANAehh8vfT5XkmSXLpS11tMVeNG+Wl3U7N69O7hdX9w05fLLLw8WNYeSlZWl++67T88//7wWLFig+fPnq2PHjjrnnHM0adIkZWRktDYyAGA/W10l/5P/aOg4uq/MuZe5CwQA7cQcc5yCN/XYsV3aViR17+kyElqh1UXNpEmTNGnSJE9fk5GRoZtuukk33XRTa+MBAA7BvvyUVLw10PD5AtPOEpPchgKAdmCyu0hH95U2bZAUmIJmKGqiFut0AkCcsmuXy77/RrBtvn6pTJ+BDhMBQPtiaefYQVEDAHHI7quW/4m/S3b/5IvuR8t840q3oQCgnYUt7bx6qey+andh0CoUNQAQh+xrz0jbNgcaxgSmnSUluw0FAO1t4FApJTWwvW+ftGa52zw4YhQ1ABBn7Po1sm+/EmybMy+UGZDnLhAAOGISk6S8Y4NtpqBFL4oaAIgjdl+1/I8/KFl/oCOnu8zF33IbCgAcCp2Cxv1qohdFDQDEEfvSk9LWTYGGMfJd/wOZlBS3oQDAodDFArSlUHbHdndhcMQoagAgTtjli2Xfez3YNl+7SGbIcIeJAMA9k9Ndym1YypmzNdGJogYA4oDds1v+J/7W0NGzD9POAGC/sCloXFcTlShqACAO2On/kkpLAo2ERPluuovVzgBgv7ApaCuWyNbWuguDI0JRAwAxzn4+V3beh8G2ufAqmd79HSYCgAgzeLiUmBTYrtorrVvpNg9ajKIGAGKYLdsh/7SHGzoG5Mmcc6m7QAAQgUxKSqCw2Y8paNGHogYAYpS1NnAdTWVFoCMlNTDtzJfgNhgARCCWdo5uFDUAEKPsh29JyxYF22bSzTLdejhMBACRK7So0cZ1srtK3YVBi1HUAEAMstuKZF/4T0PHsWNkTjvbXSAAiHTdj5Y65wSbNuRLIUQ+ihoAiDG2rk7+xx6Q9lUHOjIy5bvuDhlj3AYDgAhmjAk/W8MUtKhCUQMAMca+NUNavzrY9l17u0zHTg4TAUB0CF3a2S5fJOuvc5gGLUFRAwAxxOavlH39mWDbnDxR5vhTHCYCgCiSd6yUsH8xld0VUkG+2zxoNooaAIgRds9u+f/9F8nvD3R0zpG58ttuQwFAFDFp6dKAvGCbpZ2jB0UNAMQAa638U/8h7dge6PD55Pv2jwL/QAMAmi1sChrX1UQNihoAiAH2w1nSF58E2+aia2QGDnWYCACiU9hiAetWy9bf6wsRjaIGAKKc3bRe9rkpDR1DR8qcc5m7QAAQzY7uJ2VlB7atX1qxxGkcNA9FDQBEMVtdJf+jf5FqawIdmR3lu/luGR+/3gHgSBifT2bYqGDbLl/sLAuaj3/1ACCK2Wf/LW0pDLZ9N93F8s0A0FpDRwU37fLFsta6y4JmoagBgCjln/eh7Nx3gm3z9UvD54IDAI6IGTayobFju1S8xV0YNAtFDQBEIbu9SHbaPxs6+g2Wufhb7gIBQAwx2V2ko3oH20xBi3wUNQAQZWxtTeA6mqq9gY4O6YHlmxMT3QYDgBjCdTXRhaIGAKKMfelJqWBtsO277naZnO4OEwFA7AktarRyqWxdnbMsODyKGgCIInbJfNl3Xg22zennyJwwzmEiAIhRg4dLCfvPgO+tlDascZsHh0RRAwBRwm4rkv+xBxo6evaRueJmd4EAIIaZlFRpQF6wbVcsdhcGh0VRAwBRwFbtlf+ff5T27gl0pHSQ7zv3yCSnuA0GADGM62qiB0UNAEQ4a63s1L9LRRuDfb4b75QJWZkHAOC9sOtq1q2SrdrjLAsOjaIGACKcfecV2c/nBtvmnMtkRp/iMBEAxIk+A6S09MB2XZ20apnbPDgoihoAiGB2xRLZGVMbOoaO5H40ANBOjC9Byjs22Oa6mshFUQMAEcruLJb/0T9L1h/o6Jwj37fvkUlIcBsMAOKIGToquM11NZGLogYAIpCt2Sf/w/dLu8sDHYlJ8t32U5nMLLfBACDOhF1Xs6VQdmeJsyw4OIoaAIhA9plHw+6JYL51m0yfgQ4TAUB8Mt16SF1zg227YonDNDgYihoAiDD+j/4rO+ftYNtMOFe+U890FwgA4lzY2RqmoEUkihoAiCB23SrZZ/7V0NF/iMwVt7gLBAAIv1/NisWyfr+7MGgSRQ0ARAi7s1j+f94n1dYGOrKy5bv1JzKJSW6DAUC8yztWMiawXbFL2lzgNg8aoagBgAhgq/bI//ffS7t2BjoSEuT77o9lOnVxGwwAIJOeKYVc18jSzpGHogYAHLN1dfL/68/Spg3BPnP1rTKDj3EXCgAQJmwKGtfVRByKGgBwyFor+9y/pa8WBvvM1y+R7/SvO0wFADhQ2GIBa5bJ1uxzlgWNUdQAgEP2vddlZ89s6Dj+ZJlLr3cXCADQtP55UnJKYHvfPmntCrd5EIaiBgAcsYvnyT7/WENH30Hy3XS3jI9fzQAQaUxSkjR4eLDNdTWRhX85AcABW5Av/7//Ilkb6OicI98dv5BJSXEbDABwUOHX1XATzkhCUQMA7ayueFtgpbN91YGODmny/eBXMh07uQ0GADiksOtqNubL7i53lgXhKGoAoB3591aq9N4fNSzd7PMFlm7u2cdtMADA4R3VW+rYObBtreyKL93mQRBFDQC0E1tXp11/+ZVqN6wN9plrbpU55jiHqQAAzWWMkRk6sqGD62oiBkUNALQD6/fLTv2bqhd+EuwzZ18i3+nnOEwFAGixA+5XY+uvjYRTFDUA0MastbLP/lv209kNncedJHMZSzcDQLQJO1OzY7u0fYu7MAiiqAGANmZfeVp29pvBdvLw4+W75Ycs3QwAUchkd5ZCroO0yxe7C4Mg/kUFgDbk/+9LsjOfD7aTBg1V9s//JJPM0s0AEK3M0FHBbYqayEBRAwBtxP/RLNkZTzR09OyjTr98UL4O6c4yAQBazwwLmYK2eqmsv85dGEiiqAGANuGf96HstIcbOnK6yzf5t/JldXQXCgDgjUHHSAkJge09ldLGdW7zgKIGALxmlyyQ/c9fpfoVcbK7yHf37wPzsAEAUc+kdpD6DQ62uV+NexQ1AOAhu/JL+R+5X6rbPxUhIytQ0HTNdRsMAOApk3dscNuupKhxjaIGADxi162S/x/3SrU1gY4OafJN/q1Mj6PdBgMAeM7khVxXs3aZbE2NuzCgqAEAL9jVy+R/8FdS9d5AR3KyfN//lUyfAW6DAQDaRv8hUnJyYHvfPmndKrd54lyiFztZt26dvvzyS61du1Zr1qxRaWmpkpKS9PTTT7d4X7fffruKi4sP+viDDz6onj17tiYuAHjKLlsk/z/vDfyjJkkJifJ972cyg4a5DQYAaDMmKUkaOEzav6SzXfmlzJDhbkPFMU+KmhkzZujzzz/3YldB48ePb7I/LS3N0/cBgNawiz6T/9E/SbW1gY6kZPm+91OZ4ce7DQYAaHMmb2TwPjV25RLpoqvdBopjnhQ1gwcPVt++fTVgwAANGDBA3/nOd1q9z9tvv92DZADQdvzzPpR9/EHJ7w90pKTK9/1fygwZ4TYYAKBdmLxjZesb61fLVu0NrIyGdudJUXPxxRd7sRsAiBr+OW/LPvV/Dcs2p6XL94NfywzIcxsMANB++vSX0tID96qpq5PWLJdGjHadKi6xUAAAtJD/3Vdln/xHQ0GT2VG+H95LQQMAccb4EqTBDWfn7colDtPEN0/O1LSF1157TVu3blVSUpJ69eqlsWPHKisry3UsAHHMWis78wXZV6Y1dGZ3DtyHpkcvd8EAAM6YvGNlF38mifvVuGSstfbwT2uZSZMmeb76WUpKim688UZNnDixRfu7++67G/UlJyfr/vvvlyTVOF5TPDExUFfW1l9kjFZhPL3HmAZYa7X7qYdV+dJTwb6Ebj3U6Xd/V2L35q/IyHh6jzH1FuPpPcbUW5E2nrWF61Xy/f0LBBijblPfki+ro9tQLRQpY5qUlHTEr424MzWjR4/W8OHD1b9/f2VlZWnbtm2aPXu2Zs6cqUceeUQZGRkaO3as65gA4oitrtKuf/xRVXPeCfYl9Oytzr/9uxK6dnOYDADgWsLRfeXr1EX+0h2Stdr31RdKPeUM17HiTsSdqTmYd999V48++qh69Oihhx56yLP9btmyRW0wBM2Wk5MjSYe8Nw+aj/H0XryPqd1VKv//3SutX93QeXRf+e76nUxWdov3F+/j2RYYU28xnt5jTL0ViePp//f/k53/oSTJTDhXvmu+5zZQC0XCmBpj1KNHjyN+fdQsFDBx4kR17NhRW7Zs0fbt213HARAHbOF6+f/4w/CCZuhI+X70xyMqaAAAMSovZLGAFVxX40LUFDU+n0+5ubmSpNLSUsdpAMQ6u2S+/P/7Y2lnSbDPTDg3sGxzeoa7YACAiGOGjmxobNssG/JvB9pH1BQ1klRZWSlJSk1NdZwEQKyy1sr/35cDU86qqwKdxidz5Xdkrr5VJjHiLkUEADhmuuZKOd2DbVZBa39RU9QUFhaqqKhIKSkp6tmz+SsNAUBz2doa2Sf/ITvjPw33oEntIN/3fynfmRfIGOM2IAAgYpm8YxsaFDXtzklRM2vWLE2ePFnTp08P61+8eLHWrVvX6PkFBQV64IEHZK3VxIkTg8vOAYBXbEW5/A/+WnZuwwpn6tJNvp/8WYa7QwMADiekqLErv3S6EFU88qQ6+OKLL/Tiiy+G9dXW1urnP/95sH3ZZZfp+OOPlySVl5erqKio0bUxq1ev1owZM5STk6Pc3FxlZWVp+/btWr9+verq6jRs2DBdffXVXkQGgCC7aqn8Ux6QynY0dA4cKt9tP5PJjK57DQAA3DB5IxQsY0pLpG1FUgvuY4bW8aSoKS8v15o1a8L6rLVhfeXl5Yfdz6hRo7Rjxw7l5+eroKBAe/bsUYcOHZSXl6dx48bpjDPOkM8XNTPmAEQ4W1cn+/ozsjNfaJhuJsmcfIbMtXfItOImYACA+GKyOkk9+0ibCyQFztYYipp240lRM2HCBE2YMKHZz580aZImTZrUqH/w4MEaPHiwF5EA4JBsyTb5p/w/KX9lQ2dCosxl18ucdSHXzwAAWszkHSsbLGqWSBPOdZwofnBxCoC4Yz+fK/+T/yftrWzo7HaUfN+5R6bPAHfBAABRzeQdK/ve64HGqqWyfr8Ms4zaBUUNgLhhq6tkn5siO+ftsH5z8kSZq78rk9rBUTIAQEwYPFwyPsn6pd0V0qYNUu/+rlPFBYoaAHHBFq6X/9E/S1s3NXSmdpC55nvynTTBWS4AQOwwaelS34HS+tWSAlPQDEVNu6CoARDTbHWV7OvPyr77qlRX1/BA30HyfftHMt16uAsHAIg5Ju9Y2WBRs1Q6+xLHieIDRQ2AmGUXz5P/mUelncVh/ebrl8pcfI1MIqubAQC8ZYaOlH1rRqCxeplsba0M91hsc4wwgJhjd2wPFDNL5oc/0KWbfNfdLjPsODfBAACxb0CelJgk1dZI1XulDaulgcNcp4p5FDUAYoatrZV951XZN56V9lU3PJCQKHP2xTLnXyGTkuIuIAAg5pnklEBhs2qppP33q6GoaXMUNQBigl29TP5p/5S2FIY/MHi4fN/6nkyPXm6CAQDijsk7Vra+qFnxpXTBlY4TxT6KGgBRzW7Ml/+1ZxpPNcvsKHP5TTInTeBGmgCAdmWGjpR99elAY91K2epqZgq0MYoaAFHJblwn/+vPSIvnhT9gjMxpX5e59DqZ9Aw34QAA8a3vICm1g1S1V6qtlfJXSMNGuU4V0yhqAEQVW7g+UMws+qzxg30Gynf1d2X6D2n/YAAA7GcSEqRBx0hLP5ck2VVLZShq2hRFDYCoYDetl//1Z6UvPm38YO/+8n3jKmnkWKaaAQAighkyXDakqEHboqgBELFsXZ20dIH8H70d/LYrTK9+8l14lTTyRIoZAEBEMUNGyNY3NqyRrdork9rBZaSYRlEDIOLY4q2yc9+R/fg9adfOxk84um/gzMyoE2V8vvYPCADA4fTuL3VIk/bukerqpLUrpOHHu04VsyhqAEQEW1sju2ie7Jz/SiuWNP2knn0CZ2ZGnUQxAwCIaMa3/7qaLxdIkuzqpTIUNW2GogaAM3b/ijB28XzZz2ZLu8sbPykhUeb4k2VOO1saMoJiBgAQNcyQEbL1Rc1KrqtpSxQ1ANqVLS+VXfqF7NIF0vLFgdPyTeneU+a0s2VOniiT2bFdMwIA4AWTF3JdTcFa2ao9MqlpLiPFLIoaxCVrrbSnUtqxXdpZLLuzWNpR3LBdP/+1rjb8T39d4E/rl5JTAmvQp3QI/Ln/P1O/nZ4pZWXLZGVLmdlSVrbUsVPgOXF0Ubutq5M2rpNd+nlgFZgNaw7+5KRkmdGnBs7KDBoWV+MEAIhBR/eV0tIDxxx+v7RmhTRitOtUMYmiBnHB7iwJLKe4aqnshjVSyXapem/rdlpbG/gldeB7HaatpORAgZOVLXXqIpPdRcruInXqLNOpa2A7u0tU3nnY1tZKRRtlN+ZLBfmBPzetl/btO/iLEpOkvBEyx46VGXs6N8wEAMQM40uQBg8P3ijarvpShqKmTVDUICbZ0h2BImb1V7Irv5SKt7qO1KBmX+AM0Y7t0vrwoiesAErLkLI7Sx07yXTsHNxWx84yHTs1tJNT2vWMhvX7pYpd0s4SqbREtrREKiqULVgrbd4QKPYOp1NXmREnyBx7gpR3rExKapvnBgDABTNkhGx9UcN1NW2GogYxw27eqPJXnlL1os/kLyps3ouSk6XOOVLnbjJdcqTOXQPbGZlSQqKUkBD40+eTEkPaMtK+Kqlqr1S1V3b/n6raGzgDtHevtLtctrxMKi+TykulivLAtLXm2rM78F/RxkOf/UlMCkx1S8/Y/19W4GxHfV9ySiB7YpKUmCiz/08lJgU+izGBQmtftey+fVJNdeDMyv6+cp9P/tIS1W3ZLJWWSGU7A9PxWsLnk/rnyRx7gsyIE6SefZhaBgCIC2H3q9m4TnZPpUxaustIMYmiBlHN1tVJi+fJP/tNadVSHeSS84Cc7jJDRkiDh8sc1Uvq3E3KyPTk4Lo5e7D+Oml3RbDIsWWlgXuwlJbIlu6UynYE/ttV1rLip7YmsJ+Q+7k0mvIWmqP5e5akQ49pUxISpZ69ZfoMlHr3D/zZs49McvRNpwMAoNV69pEyMgPHANYvrVkujRzjOlXMoahBVLLlZbJz3pb9cFbg7EFTuubKDBkuDTlWZvDwwJkYh4wvoeFaGvU9aCFk6+oChU/ZDqlsp+yundKu0v3b+wuhsp2BKWC2pSWKhzqkS526BK4L6por9R4QKGCO6i2TlOQuFwAAEcT4fIHrar74VNL+62ooajxHUYOoYa2V1q+Wnf2m7Odzm7x2I7HfIKWdc4l29xks06Wbg5StZxISgsWCdPCzQLauLnBfl8oKqXK3VFkhW7lbqiwPtlW5W7ZmX+BsTm1t+J91tYFtawOLFySnBKbj7d82ySlSUrJSszoqoVNnVSZ3kOncVerUNVDIsCQlAADNYoaMkA0WNV85ThObKGoQFWzxVvmn/0v6amHjBxMSAssAn3G+upx0mowxqiwubv+Q7cwkJOxfOKBTQ18bvE/HnMAZrr1xMKYAALSFsOtqCtfJVu5mtU+PUdQgotm6Otl3X5N97enGywJ37Cwz/pzADRqzO0sSF58DAIDIc1RvKbNjw9TxNV9Jo05ynSqmUNQgYtkNa+R/8h9S4frwBwbkyZx5ocxxJ8kk8iMMAAAimzFGZvBw2YUfSwpMQTMUNZ7iiBARx1btlX31adn33ghfBSyzo8ykm2VOHM8ZGQAAEF2GjJDqixruV+M5ihpEFLtkvvzTHwnc2DGEOfUsmW/eIJOR5SgZAADAkTN5IdfVbFovu7uc4xoPUdQgItiqPbJT/xFY1SxUbk/5rr0tcH8ZAACAaNX96MBtHcrLAu3Vy6TjT3aZKKb4XAcAbNlO+f/8s/CCJiFR5oIr5Pv1QxQ0AAAg6hljwo5p7CqmoHmJMzVwym4plP+h30o7tjd0DhwWODtzVG93wQAAALw2ZIS0YI4kihqvUdTAGbt6mfz/d6+0Z3ewz5x9scxlNwTuvgsAABBDwu5Xs7lAtmKXTGZHl5FiBkeOcMJ+Plf+B3/VUNAYI3PFLfJdfhMFDQAAiE25R0kdOze0V3/lLkuM4egR7c7/zqvyP/pnqbYm0JGYJN93/0e+sy50GwwAAKANNbquhqWdPUNRg3Zj/X75n5si+/xjgbvpSlJahnx3/15m9KluwwEAALSHPBYLaAtcU4N2YWv2yf/YA9LCTxo6u3ST787fyPQ42l0wAACAdmSGDG+4rmZLoWx5qUxWJ5eRYgJnatDmbF2d/P+8L7yg6T1Avp/+mYIGAADEl5weUqeuwaZdxXU1XqCoQZuzz02RvlrY0HHMcfLdc69MR76VAAAA8eXA62rEFDRPUNSgTfnff0N29psNHaNOku+OX8qkprkLBQAA4NKQ4cFNrqvxBkUN2oz9aqHss1MaOvoMlO+WH8okcikXAACIX2FnarZuli3b4S5MjKCoQZuwmwvk/9efJOsPdHTqKt8dP5dJSXEbDAAAwDGT013q0i3Y5rqa1qOogedsean8f/+9VLU30JGSKt8dv5DJ7uI2GAAAQIQwgxumoHFdTetR1MBTtmaf/P/3R2nH9kCHMfJ9+0cyvfu7DQYAABBJQu9Xs3qZwyCxgaIGnrHWyj7xN2ndqmCf+eaNMiPHOkwFAAAQecLO1GzbLFu2012YGEBRA8/Y15+Rnf9RsG1O/7rM1y5ymAgAACAyma65UuecYNuu5rqa1qCogSf88z6Uff3Zho6hI2Wu+q6MMe5CAQAARDATsrSzKGpahaIGrWY3bQhMO6vX/Wj5bv0xSzcDAAAcyuDQ+9VQ1LQGRQ1axdbVyf/E36TamkBHRqZ83/+lTFqG22AAAAARLuy6mq2bZMtL3YWJchQ1aBX7zitSwdpg23fz3TLdergLBAAAEC1yukuduja0WQXtiFHU4IjZrZtkX50ebJtTzpQZPtphIgAAgOhhjJEZfEywzRS0I0dRgyNi/X75p/69YdpZx04yk252GwoAACDahF5Xw2IBR4yiBkfEzp4prV0RbPuu+Z5MOtfRAAAAtIQZ0nATThVtlK3Y5S5MFKOoQYvZ4q2yLz8ZbJsxp8kcd5LDRAAAAFGqWw+pY+eGNtfVHBGKGrSItVb+p/5Pqq4KdGRkylz1HbehAAAAolSj62qYgnZEKGrQInbuO9KKJcG2ufI7MpkdHSYCAACIciFT0ChqjgxFDZrNlu6QfeHxho6RY2XGnu4uEAAAQAwIu1/Npg2yu8vdhYlSFDVoFmut/NP+Ke3dE+jokB5YHMAYt8EAAACiXfeeUlZ2Q3vNcmdRolWiFztZt26dvvzyS61du1Zr1qxRaWmpkpKS9PTTTx/R/iorK/XCCy9o/vz5KisrU3Z2tsaMGaNJkyYpPT3di8hoITv/I+nLBcG2ufxGmU5dHCYCAACIDYHraobLfj5XkmRXLWURphbypKiZMWOGPv/8cy92pYqKCv3iF7/Qli1blJubqzFjxmjTpk166623tGjRIt17773KzMz05L3QPLa8TPbZRxs6ho6UGfc1d4EAAABizeDhUn1Rw3U1LeZJUTN48GD17dtXAwYM0IABA/Sd7xz5alhTp07Vli1bNHbsWN11111KSEiQJD3++OOaNWuWpk6dqjvuuMOL2GgmO+M/0u6KQCMlVb7r7mDaGQAAgIfM4OGy9Y1NG2Qrd3MPwBbw5Jqaiy++WJMmTdLo0aOVnZ19xPspKyvTnDlzlJCQoFtuuSVY0EjStddeq6ysLM2dO1dlZWWtD41msZs3yn72QbBtLr1Opmuuu0AAAACx6KheUkZWYNtaaQ33q2mJiFooYNGiRbLWatiwYY2Ko6SkJI0ePVp+v1+LFy92ki8e+V97OvAXS5J69JKZcK7TPAAAALHIGBOYgrYfU9BaJqKKmoKCAklSv379mny8vn/Dhg3tFSmu2YK10hefBtu+i66W8SUc4hUAAAA4UiasqOFMTUt4ck2NV0pKSiRJnTt3bvLxLl26hD2vOe6+++5GfcnJybr//vslSV27dm1pTE8lJgb+F+Tk5DjN0ZSdD/9R+/ZvJ/YfrC5nXyjji6g6uJFIHs9oxZh6i/H0HmPqLcbTe4ypt2J5PGtOOk076hdnKlynLmkd5GuH62piYUwj6gi1qqpKkpSSktLk4/X91dXV7ZYpXu1bvkT7vvgs2M68+rsRX9AAAABEs8Te/WUy919X4/dr34olbgNFkYg6U2P3X7vh5cpaDzzwwCEfLykpCb6vC/UVcXFxsbMMB7LWyv/E3xs6BuRpV++BMhGU8WAicTyjHWPqLcbTe4yptxhP7zGm3or18bQDhkmLA18s71rwsSr6DG7z94yEMTXGqEePHkf8+oj66r1Dhw6SGs7YHKj+DM3BzuTAIysWSyHzOH2XXMsSzgAAAO3ADAm5rmYViwU0V0QVNfXXt+zcubPJx3fs2BH2PHjPWiv/y9MaOoaOlBkywl0gAACAOBK6WIA25stW7XEXJopEVFHTp08fSdL69eubfLy+v/55aANL5kkb1gSbvou/5TAMAABAnDm6j5SWHtj2+6W1K9zmiRIRVdSMGjVKxhitWLFCu3btCnuspqZGCxculDFGxx13nKOEsc36/fK/8nRDx8ixMv2HuAsEAAAQZ4wvQRp0TLDNFLTmcVLUzJo1S5MnT9b06dPD+jt16qRTTz1VtbW1mjJliurq6oKPTZs2TeXl5Ro3blyjG3PCG3bBHGlzQbDtu+gah2kAAADik+EmnC3myepnX3zxhV588cWwvtraWv385z8Pti+77DIdf/zxkqTy8nIVFRWptLS00b5uuOEGrVmzRvPmzdPkyZM1YMAAFRYWqrCwULm5ubr++uu9iIwD2Lo62deeCbbNmNNkejV9E1QAAAC0HTNkuIJr8xasla3aK5PawWWkiOdJUVNeXq41a9aE9Vlrw/rKy8ubta+srCzdd999ev7557VgwQLNnz9fHTt21DnnnKNJkyYpI6Ptb0AUj+wn70nbiwIN45O58Cq3gQAAAOJVr35ShzRp7x6prk7KXykdw+UXh+JJUTNhwgRNmDCh2c+fNGmSJk2adNDHMzIydNNNN+mmm27yIB0Ox9bUyL7xXLBtTjlDpvvRDhMBAADEL+NLkAYOk5Z+LikwBc1Q1BxSRC0UADfsnP9KO/ffbCkhUeaCK90GAgAAiHNh96vhuprDoqiJc3ZftezMF4Jtc9rZMl1zHSYCAACAGRxyn8D1a2T334QeTaOoiXN23ofSrv0LNiQly5x/udtAAAAAkHr3l1L2Lw5QVyutW+k2T4SjqIlj1lrZ2W8G2+bUM2WyuzhMBAAAAEkyCQnSwLxg265e5jBN5KOoiWf5K6XC9cGmmXC+uywAAAAIE3a/mjUUNYdCURPH7OyZDY0hI2R69nYXBgAAAGHM4GMaGutWydbUuAsT4Shq4pQtL5Vd+HGw7TuDszQAAAARpe8gKSk5sF2zT9qw5tDPj2MUNXHKfvR24KIzSerUVRp1ottAAAAACGMSk6T+Q4JtlnY+OIqaOGTr6mQ/nBVsm9O/HrgYDQAAABEl7LoaFgs4KIqaeLR4nlS2I7CdkChz+tlu8wAAAKBJYdfV5K+QratzFyaCUdTEIX/oMs6jT5XJ6uQwDQAAAA6q3xApITGwXV0lbcx3mydCUdTEGbt5o7RqabBtzjjPYRoAAAAciklJkfoNCraZgtY0ipo4Yz9oOEuj3v2lAXkHfzIAAACcM4MapqCxWEDTKGriiN1TKfvp7GDbTDhPxhh3gQAAAHBYoYsFaM1yWT/X1RyIoiaO2E9nB+ZiSlJahszY8W4DAQAA4PAG5km+/YfteyulzRvd5olAFDVxwlobNvXMjDsrMEcTAAAAEc2kpkm9BwTbTEFrjKImXqxYIm3dHNg2Rmb8uW7zAAAAoNlCl3ZmsYDGKGrihH/2zIbG8NEy3Xq4CwMAAIAWCV0sQGuWyVrrLkwEoqiJA3ZHsbRkfrDtO+N8h2kAAADQYoOOkeoXeKrYJW3d5DZPhKGoiQP2w7ck6w80crpLxxznNhAAAABaxKRnSD37BNt2FdfVhKKoiXG2Zp/snLeDbTPhPBkf/9sBAACiTfjSzlxXE4qj2xhnP/9Y2l0eaCQny5x6lttAAAAAOCIHLhbAdTUNKGpinP343eC2OXFC4NQlAAAAok/oYgFlO6Tire6yRBiKmhhmS3dIIeuYc5YGAAAgepmsbKn70cG2ZQpaEEVNDLML50r1pyW75kr9h7gNBAAAgFYJnYImFgsIoqiJYXb+nOC2GXOaTP0ygAAAAIhOIYsFcKamAUVNjLLFW6X1q4NtM/Y0h2kAAADghbCbcJZsk91Z7C5MBKGoiVF2QcNZGvXoJfXs6ywLAAAAvGE6dw3cd3A/u5qzNRJFTcyy8z8KbpuxTD0DAACIFWFna1ZzXY1EUROT7OaN0uaCYNuMOd1hGgAAAHiK62oaoaiJQfbzkKlnfQbK5B7lLgwAAAA8FbYC2tbNsuWl7sJECIqaGGOtDZ96NoYFAgAAAGJK11wpu0tDm+tqKGpizsZ8afuWYNOcMM5hGAAAAHjNGCMTOgWNooaiJtaEnqXRwGEyXXLchQEAAEDbCJmCZlksgKImlli/X3bB3GCbe9MAAADEptAzNdpcIFtZ4S5MBKCoiSX5K6XSksC28cmMPtVtHgAAALSN7j2lzI4N7ThfBY2iJoaETT0beqxMVrazLAAAAGg7xpjwKWirKGoQA2xdnezCj4NtVj0DAACIbWYQ96upR1ETK1Z+KVXsCmwnJMocf7LbPAAAAGhTYfer2bhOtmqPuzCOUdTECLsgZOrZ8ONl0jLchQEAAEDb69lbSksPbFu/tHal2zwOUdTEAFtTI/vFZ8G2GXu6wzQAAABoD8aXIA0cFmzH8xQ0ippYsGyhtLcysJ2cIjNyrNs8AAAAaBcm7H41FDWIYnb+nOC2GTlWJiXVYRoAAAC0FzMo5LqaDatl91W7C+MQRU2Us9VVskvmB9vccBMAACCO9B4gJacEtmtrpfVr3OZxhKImytkl86X6irxDunTMaLeBAAAA0G5MYqI0IC/Ytmu+cpjGHYqaKGcXhEw9O+4kmaQkh2kAAADQ3riuhqImqtk9ldJXC4NtVj0DAACIP6E34VT+StnaWndhHKGoiWYrlgTmTkpSRqaUd6zbPAAAAGh//QZJiYmB7X3V0sZ8t3kcoKiJYnbp58Ftc8zxMgkJDtMAAADABZOcIvUbHGzH4/1qKGqilLVWNmTqmUac4C4MAAAAnAqdghaP19VQ1ESrwnXSrtLAtjEyxxznNg8AAACcCV0sQGuWy/rr3IVxgKImStmlIWdp+g2WychyFwYAAABuDRgi+fYf2u+tlDZvdJunnVHURKnQqWdmBPemAQAAiGcmNS1wI8794m0KGkVNFLKVFVL+qmDbcD0NAABA3Au7X02c3YSToiYK2eWLJesPNDI7Sr36O80DAAAA98ygkOtqVi+TtdZdmHZGURONQpdyHj5axsf/RgAAgLg3aFjDdsUuadtmd1naGUfDUcb6/bJffdHQwdQzAAAASDLpmVLPPsF2PF1XQ1ETbQryA5W3JPl8MsNGOY0DAACAyBG+tDNFDSJU2A03++fJpGe4CwMAAIDIEqc34Uz0akf79u3TK6+8oo8//lglJSXKyMjQyJEjdcUVV6hLly7N3s/tt9+u4uLigz7+4IMPqmfPnl5Ejko29HoalnIGAABACDNomILLA+wslt2xXaZLN5eR2oUnRc2+ffv0+9//XqtWrVKnTp10wgknqLi4WB988IG++OIL/eEPf1D37t1btM/x48c32Z+WluZF5KhkK3ZJG9YE2yzlDAAAgFAmu7PU7Shpe5GkwNkaczJFTbO8/PLLWrVqlQYPHqxf/OIXSk1NlSS98cYbevLJJ/Xwww/rt7/9bYv2efvtt3sRLabYZYuk+qX5sjtLR/d1mgcAAACRxww+RnZ/UaM1y6STz3AbqB20+pqa2tpazZo1S5J08803BwsaSbrgggvUp08frVixQuvWrWvtW2Fpw/U0ZvhoGWMchgEAAEBECrlfTbxcV9PqomblypWqrKxUbm6u+vXr1+jxE088UZL0+eefN3oMzWf9dbLLGpZy5noaAAAANCVsBbRtm2V3lboL005aPf2soKBAkposaCSpf//+Yc9rrtdee01bt25VUlKSevXqpbFjxyorK6t1YaPZ+jVSZUVgOyFBGjrKaRwAAABEqC7dpM5dpZ0lgfaaZdIJ49xmamOtLmpKSgKDdbAVzjp37hz2vOaaNm1aWHvq1Km68cYbNXHixBbt5+67727Ul5ycrPvvv1+S1LVr1xbtz2uJiYH/BTk5OYd8XsU7L6ty/3by0JHq3LvPIZ8fr5o7nmg+xtRbjKf3GFNvMZ7eY0y9xXg2T9mI0ar68L+SpNTCdco695KDPjcWxrTVRU1VVZUkKSUlpcnH66+xqX/e4YwePVrDhw9X//79lZWVpW3btmn27NmaOXOmHnnkEWVkZGjs2LGtjR11qhd+GtxOGX2ywyQAAACIdMnDRgWLmn3LF7sN0w5aXdRYa1v1+IFuuummsHavXr103XXX6aijjtKjjz6qp59+ukVFzQMPPHDIx0tKSlqc0Uv1FfGh7s1jd5XKn78y2K7sN1R7DvH8eNac8UTLMKbeYjy9x5h6i/H0HmPqLcazeWyPhlk9tQX52r5hnUx6ZpPPjYQxNcaoR48eR/z6Vi8U0KFDB0lSdXV1k4/X94euinYkJk6cqI4dO2rLli3avn17q/YVbUIXCFDnrtJRvdyFAQAAQOTr3lPK7BjYtlZau8JtnjbW6qKm/pqUHTt2NPn4zp07w553pHw+n3JzcyVJpaWxv4JDmLClnE9gKWcAAAAckjEmrpZ2bnVR06dP4NTW+vXrm3y8/v409c9rjcrKwKXyrT3rE01sXV3gppv7sZQzAAAAmsMMHh7ctmsoag4pLy9PaWlp2rZtW5OFzbx58yRJxx9/fKvep7CwUEVFRUpJSVHPnj1bta+okr9S2rt/3bPERCnvWLd5AAAAEBXC7ldTsFa2aq+7MG2s1UVNYmKizjnnHEnS448/HrbK2RtvvKGCggLl5eVp4MCBwf5Zs2Zp8uTJmj59eti+Fi9eHDyzE6qgoEAPPPCArLWaOHFicNm5eGC/aph6pkHHyKR2cBcGAAAA0aNnbyktPbDt9we+LI9RnlQHl156qZYuXapVq1bpzjvvVF5enkpKSrRmzRplZmbqtttuC3t+eXm5ioqKGl0bs3r1as2YMUM5OTnKzc1VVlaWtm/frvXr16uurk7Dhg3T1Vdf7UXkqGFDr6cZcYLDJAAAAIgmxpcgDRwmfblAUuC6GnPMcY5TtQ1Piprk5GT9+te/1ssvv6y5c+dqwYIFSk9P1/jx43XFFVc0e5GAUaNGaceOHcrPz1dBQYH27NmjDh06KC8vT+PGjdMZZ5whn6/VJ5eihi3dIW1qmNJnhnM9DQAAAJrPDD5Gtr6oWfOV4zRtx7N5XMnJybriiit0xRVXHPa5kyZN0qRJkxr1Dx48WIMHD/YqUtQLm3rWNTewNB8AAADQTGbQMQrekXH9atmafTJJyS4jtYn4Oe0RjVYsCW6a4aNZyhkAAAAt03uAlJwS2K6tldavdpunjVDURChrbdh64mYoq54BAACgZUxiojQgL9iO1fvVUNREquIt0q6dDe2QmycBAAAAzRW6tHOs3q+GoiZChVXRPXrJZHZ0FwYAAABRywxquAmn8lfK1ta6C9NGKGoiVejUM87SAAAA4Ej1GxS4ibskVVdJG/Pd5mkDFDURKuzU4GCKGgAAABwZk5wi9WtYYTgWp6BR1EQgu7NYKtkWbHOmBgAAAK0ROgUtFhcLoKiJQHbN8oZGTneZzs27eSkAAADQlNDFArRmuay/zl2YNkBRE4lWN9ztlbM0AAAAaLUBQyTf/kP/vZXS5o1u83iMoiYChZ0S5HoaAAAAtJJJTQvciHO/WJuCRlETYWx5mbR1U7DNmRoAAAB4Ifx+NV8d4pnRh6Im0oReT5PdWcrp7i4LAAAAYkbYl+Wrl8la6y6MxyhqIkzoEntm8HAZYxymAQAAQMwYNKxhu2KXtHWzuyweo6iJMDZkkQAx9QwAAAAeMemZUs8+wXYsTUGjqIkgtnK3tGlDsG1YJAAAAAAeMoMb7lejGFosgKImkqxdIdXPbczIknr0cpsHAAAAMSVssYAYuq6GoiaChJ0CHDSM62kAAADgrdDLG0pLpB3b3WXxEEVNBAldL5ypZwAAAPCa6dhJyu0ZbMfK/WooaiKErdorbcwPtsPmOwIAAAAeCfvyfA1FDby0bpVUVxfY7pAmHd3XaRwAAADEqEGh19XExgpoFDURIvT+NBo4TMaX4C4MAAAAYlbYmZrtW1S3s9hdGI9Q1ESI0CrZcH8aAAAAtBHTpZvUOSfY3rd8icM03qCoiQB2X7W0bnWwzSIBAAAAaEuhx5s1yxY5TOINipoIULNmhVRbE2gkJ0t9BrgNBAAAgNgWsijVvmWL3eXwCEVNBNi3fHFDY8BQmcQkZ1kAAAAQ+0Ivd6jduE7+8l0O07QeRU0ECK2OuZ4GAAAAbS73KCkrO9jctyK6r6uhqHHM1tWqZuWXwTbX0wAAAKCtGWPCvkyP9iloFDWO1a5bHbjxpiQlJEr9BrsNBAAAgPgwOLSoie7FAhJdB4h3YT9A/QbJJKe4CwMAAICIY62V30p2/3bgT+3v2/9Y6ONW8oe8TpL8+/ut9r/OWvl7DlVdWq6sMfJvr1CX0gp16pTp7HO2BkWNY2HX04SsQgEAAGLXgQepDQeeBxykhh7Ahr5u//MOPEhtONhtfODb8Nz926Ht0NcFc0l+NTzXr8AD/gP3o4b3rs/mtzbs8zTKEfJ56rcV8n4Hfp7Q90vpUC5rrfbs2dvk5w0fv5DP0MS4NzrID3m/YLaQfeqAsfbbJnI39X6hn7f+PULeL7Dv8P+3oYVKmxr7w+DmjYvX6eIzRrb1O7YJihqHrN8fdrMjFgkAgOjV1EFq0wdZLTtIDR7sNfMgNfwAruUHqRklflkrlZeXNzpIDf08jQ72WnGQGvq6Aw9SGx3sHeYg9cADzyM5SG3qm/CmDjwb5276IFVaJf8B/y/QGqWuA8Qs//ZtriMcMYoal4oKZCsrAtvGJw3Mc5sHQJiWHKQ2HMAd+iC1MmGPrLXaUVZ90AOn1h6k2uB7HzjdIHAo1ZqD1EaZww48m3+QGv5t5REepO7fR2JikfxW2lezr/EUjNDXhe7nCA9S7QFj1m7fpLarLa4DAGgjRpIxks8EWj4jmbo6mdoa+YyUkJ7sOOGRo6hxyK5e1tDo3V8mNc1dmDhw4DeEoQdOBx7sNfXt3uG+SbUHHPA0/iaweQepYQeaHh+kNvq2MuSzSU0fpKakVshaq8q9e4/oIDXsoLHRWDYev4NOwVB4tkMdpDb1TXj9tuLyIDXW7XUdAHCmyYNUBVa28pnAY8H2/ufWb8tof1/oc03wNT7T8Jjq2/sf9x3kuWb/PhXyfsHH1PAe9e/nC9mHCXm/pp5bnzstrYOMMarauzfs85oD3s+nAz5TSO6GzxT6GcKz+MIyh4/Tgbl9Ye9vGsYr9LlNjN+B/4/qx6/R5wl5XWhuX8jrDszSeNwD7QPZ3eXq0qmTErI7q7i4uNU/k65Q1DhS67f6eF2ZbLdR8stIA0dL63Y1/gY27CA55AD2CA9SD3qwF3rgHfrc0CyeH6Q2URyEvN+BB6kHfhPe1EGqtEZ+a0PGkYNUANGlqQOZxgdUCjtoPJKD1NCDxtD3S0lOkoxUV1Nz6IOsIzhIPfBgL/zA88BsBxxMe3yQeugDz8bPDdtPCw9Su3TuLJ+RykpLG94z9HWh7xXyurDPpOYfpMa6nJwcSYrqA/BIYjKylJDd2XWMVqOocaSmzq8HOpwgDTsh0FEn6VNO+QOxKPRAMHAgZWRkm3WQ2tQBbfgBXOMDwQMPig987uEOUsO+8W3iYO/Ag9RGB7BN5Q7ZDssd8rqmcjc+KG98sJeVlSmfMdpdURGy7yM/SA072D3MN6kHLQ6i+CCVA0bv5XQOzMToUFvpOAkQuyhqHDHbi1xHQCscePB1sG9Sww5WQg7gDjxoPJKD1EMdZIXup6kDp5YcpKalpcknqapq7wEHewceaDc+SA1+c9zCg9SwA8/g/g448Dzg2+D2OEg98Jvwpn4Gmhr3UBwweq9hTN0WAwAAdyhqHPElJmqQyuWr2isjK1+vfo2/rdSRHaQeOK/2UN+kHmq+qBcHqY0PykO+CW7lQWrogWd9ji6dO8sYqbS0tNkHqY0Oblt4kBrrOAgHAACRjqLGkeTcHvrLNT2Uk5MjW7NPJWW7XEeKCTld0yVJxf49jpMAAACgvfhcB4BkkqJ3+TwAAADANYoaAAAAAFGNogYAAABAVKOoAQAAABDVKGoAAAAARDWKGgAAAABRjaIGAAAAQFSjqAEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENUoagAAAABENYoaAAAAAFGNogYAAABAVKOoAQAAABDVKGoAAAAARDWKGgAAAABRjaIGAAAAQFSjqAEAAAAQ1RJdB3DNGOM6gqTIyRErGE/vMabeYjy9x5h6i/H0HmPqLcbTey7HtLXvbay11qMsAAAAANDumH4GAAAAIKpR1Dj2k5/8RD/5yU9cx4gZjKf3GFNvMZ7eY0y9xXh6jzH1FuPpvVgY07i/psa1ffv2uY4QUxhP7zGm3mI8vceYeovx9B5j6i3G03uxMKacqQEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBU4z41AAAAAKIaZ2oAAAAARDWKGgAAAABRjaIGAAAAQFSjqAEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENUSXQeIJfv27dMrr7yijz/+WCUlJcrIyNDIkSN1xRVXqEuXLi3aV2VlpV544QXNnz9fZWVlys7O1pgxYzRp0iSlp6e30SeIPF6N6fLly7Vs2TKtXbtWa9euVUVFhY466ij99a9/bbvwEciL8aysrNSiRYu0cOFCbdiwQSUlJTLG6Oijj9a4ceN09tlnKzExfn61eDGmdXV1evHFF5Wfn6/NmzervLxcdXV16tKli4499lhdfPHF6tq1axt/ksjh5e/SUFu2bNGPfvQj1dTUaOTIkfr5z3/uYerI5dV43n777SouLj7o4w8++KB69uzpReSI5/XP6NatW/XKK69o6dKlKisrU2pqqnr06KGxY8fqwgsvbINPEFm8GM8PPvhA//znPw/7vNtvv13jx49vbeSI5+XP6OLFizVz5kzl5+drz549Sk9P18CBA3X++edrxIgRbfQJWs5Ya63rELFg3759+v3vf69Vq1apU6dOysvLU3FxsdauXausrCz94Q9/UPfu3Zu1r4qKCv3iF7/Qli1blJubq/79+2vTpk0qLCxU9+7dde+99yozM7ONP5F7Xo7pPffco4KCgrC+eCtqvBrPZ599Vi+99JKMMerXr5+6d++u8vJyrVq1SjU1NcrLy9PPf/5zpaSktMOncsurMa2qqtJ1112n1NRU9enTR506dVJtbW2waExLS9OvfvUr9e/fvx0+lVte/r0/0G9/+1stX75c1tq4KWq8HM/6ouZgB4RXX321OnXq5GX8iOT1z+j8+fP10EMPqba2Vn379lWPHj20e/dubdy4USkpKfr73//ehp/GPa/Gc+XKlXrvvfeafGzPnj1asGCBJOnvf/+7cnNzPf0MkcbLn9E33nhDTz75pIwxGjJkiDp37qxt27YpPz9fknTLLbfo7LPPbsuP02zx83VqG3v55Ze1atUqDR48WL/4xS+UmpoqqeGH4eGHH9Zvf/vbZu1r6tSp2rJli8aOHau77rpLCQkJkqTHH39cs2bN0tSpU3XHHXe02WeJFF6O6ciRI3XyySdr4MCByszM1I9//OO2jB6RvBrP1NRUXXLJJfr617+uzp07B/u3bNmi3//+91q5cqVefPFFXX311W32WSKFV2OalJSk3/3udxo0aFDw77sk+f1+Pfvss3rllVf02GOP6d57722zzxIpvPx7H+r999/XsmXLdNZZZ+ndd9/1OnbEaovxvP3229siatTwckw3bNigv/71r+rQoYPuuece5eXlBR/z+/1av359m3yGSOLVeObl5YWNX6i3335bCxYs0JAhQ2K+oJG8G9Py8nJNnz5diYmJ+tWvfhU2vp999pkefPBBPfXUUzr99NOD7+ES19R4oLa2VrNmzZIk3XzzzWH/Yy+44AL16dNHK1as0Lp16w67r7KyMs2ZM0cJCQm65ZZbwg5wrr32WmVlZWnu3LkqKyvz/HNEEi/HVJK+9a1v6dJLL9Wxxx4bV9P36nk5nhdffLGuuuqqsIJGknr06BEsZD7++GMP00cmL8c0ISFBeXl5YX/fJcnn8+mKK65QUlKS1qxZo6qqKm8/RITx+u99vV27dumpp57SiBEjdOqpp3qaOZK11XjGM6/H9D//+Y9qa2t12223NTog9/l8GjBggHfhI1B7/YzOmTNHknT66ae3aj/RwMsxXbNmjWprazV8+PBGP58nnXSSevfurerqam3atMnbD3GEKGo8sHLlSlVWVio3N1f9+vVr9PiJJ54oSfr8888Pu69FixbJWqthw4YpOzs77LGkpCSNHj1afr9fixcv9iJ6xPJyTNF+49m3b19JUmlpaav2Ew3aa0yNMfL5fDLGNCp6Yk1bjel//vMf7du3T9/+9rc9yRkt+D3qPS/HdNOmTVqxYoV69Oih0aNHe541GrTHz+j27du1atUqJSYm6uSTTz7i/UQLL8c0KSmpWe+ZkZHRspBthOlnHqi/VqOpHx5JwXnwB17TcST76tevn2bPnq0NGzYcQdLo4eWYov3Gc9u2bZLUqCCPRe0xptZavfLKK6qurtaIESOa/Q9MtGqLMf3iiy/0ySefaNKkSerevbt27NjR+qBRoq1+Rl977TVt3bpVSUlJ6tWrl8aOHausrKzWhY0SXo7pV199JUk69thjtW/fPn3yySfBb8/79Omjk08+WWlpaV7Ejljt8Xv0o48+kiQdf/zxEXPw3Za8HNMBAwYoLS1NX331lVauXBl2tmbevHnauHGjhgwZcsTXOXqNosYDJSUlknTQ1STqp+nUP685+zpwak+9+vdozr6imZdjivYbz5kzZ0qSTjjhhFbtJxq01ZhOmzZNu3bt0t69e1VQUKBt27apZ8+e+u53v9u6wFHA6zGtqqrSY489pqOOOkoXX3yxJxmjSVv+jIaaOnWqbrzxRk2cOPEIUkYXL8e0sLBQkpScnKz/+Z//UVFRUdjj06dP1w9/+EMNGzasNZEjWnv82zR37lxJ8TH1TPJ2TNPT03Xrrbfqb3/7m379618HFwrYvn278vPzNWrUKN12223ehW8lihoP1M9zP9hqT/XzGZszH/5w+6rvr66ubnHOaOLlmKJ9xvPtt9/W0qVLlZ6eHhcHkG01pvPmzQue8ZKkXr166Qc/+IG6det2hEmjh9dj+uyzz6q4uFi/+tWv4mqZ8Xpej+fo0aM1fPhw9e/fX1lZWdq2bZtmz56tmTNn6pFHHlFGRobGjh3rTfgI5eWYVlZWSgp8GZSenq4f/ehHGj58uMrKyjRjxgzNnTtXf/7zn/XAAw/E7Kpybf1v09q1a1VUVKT09HQdf/zxRxYyyng9pieddJIyMjL04IMPauXKlcH+jh076phjjomo1Xjj77d8GzjcqtgtWTW7/rnGmFZlinZejinafjyXL1+uJ554QsYYfe973zvomcZY0lZjWr98a3l5udatW6dnn31WP/nJT/Td735XEyZMOKJ9RgsvxzQ/P1+zZs3S6aefruHDh7c2WlTy+mf0pptuCmv36tVL1113nY466ig9+uijevrpp2O+qPFyTP1+v6TAfaq+//3va+TIkZKktLQ0/eAHP9CWLVuUn5+v//73v7ryyiuPPHQEa+t/m+qnnp1yyilx88WG12P6+uuva9q0acF7JXbr1k3bt2/Xc889p2nTpmnNmjX64Q9/2JrInmGhAA906NBB0sHPntT3N2e5u/p9HayCrt9XrN8DxMsxRduOZ0FBgf785z+rtrZWN9xwQ8wf1NRr65/RrKwsjRo1Sr/61a/UqVMnTZkyJeanW3o1pnV1dfrXv/6ltLQ0XXfddd6GjCLt9Xt04sSJ6tixo7Zs2aLt27e3al+RzssxrX9O586dgwVNqDPOOEOStGzZsiPKGg3a8me0rq5On376qaT4mXomeTumy5cv11NPPaW+ffvq7rvvVu/evZWamqrevXvrhz/8ofr166d58+ZpyZIl3n2AVoiPsrWN1d/p+2AXoO7cuTPsec3ZV/1rDlT/HrF+d3EvxxRtN55bt27Vvffeq8rKSl1++eU699xzWxc0irTXz2haWpqOP/54vf322/ryyy9j+roFr8Z0x44d2rBhg7Kzs/XAAw+EPVY/5Wft2rX6zW9+o9TUVP3kJz9pbfSI1F4/oz6fT7m5udq1a5dKS0tjeqqkl2NaP045OTlNPl7fX15e3uKc0aItf0aXLFmiXbt2KTc3V0OGDDnykFHGyzH98MMPJQVWTPP5ws+D+Hw+jR07VuvXr9eyZcuaLMzbG0WNB/r06SNJB71JVuhqJq3dV31/c/YVzbwcU7TNeO7cuVN/+MMfVFZWpvPOO0+XX35564NGkfb8Ga1fWSqWD24k78e0rKzsoPf0qqys1PLly2N6dan2/BmtLxZj/ey5l2NavwT+7t27m3y8oqJCUmyPaVv+jNbfm+a00047wnTRycsxrS+A6s/+HKi+/2A/w+2NosYDeXl5SktL07Zt27R+/fpGy+jNmzdPkpp1kdqoUaNkjNGKFSu0a9cudezYMfhYTU2NFi5cKGOMjjvuOG8/RITxckzh/Xju3r1b9957r7Zv364JEybo+uuv9zxzpGvPn9Hly5dLUsQsm9lWvBrTbt266fnnn2/ysWXLlum3v/2tRo4cqZ///OfeBI9Q7fUzWlhYqKKiIqWkpKhnz56t2lek83JMR4wYoZSUFG3dulUlJSWNvjmv/3t/sKV5Y0Fb/YxWVVUF78MSb0WNl2Nafwyan5/f5OP1/ZFydpZrajyQmJioc845R5L0+OOPh10P88Ybb6igoEB5eXkaOHBgsH/WrFmaPHmypk+fHravTp066dRTT1Vtba2mTJmiurq64GPTpk1TeXm5xo0bF/P3AfFyTOHteFZXV+u+++5TYWGhTj75ZN16661xubCFl2P6+eefB2+8G6q6ulrPPPOMli9fruzsbI0aNartPlAE4O+9t7wcz8WLFzd5B/KCggI98MADstZq4sSJMX8xtpdjmpKSonPPPVd1dXWaMmVK2L4WL16sDz/8UMYYnXXWWW38qdxpq7/z8+bNU3V1tQYNGqQePXq03QeIQF6Oaf01snPnzm10s84FCxZo7ty5MsZEzLW0sf3bpx1deumlWrp0qVatWqU777xTeXl5Kikp0Zo1a5SZmdloHe/y8nIVFRU1eef1G264QWvWrNG8efM0efJkDRgwQIWFhSosLFRubm7cfCvu5Zi+9957ev/99yUFznhJUnFxcdg3tTfffHPwplSxyKvxfOaZZ7RmzRr5fD4lJCTo4YcfbvL9br/99jb7LJHCqzFdt26dZsyYoU6dOqlfv35KS0tTWVmZNmzYoN27dystLU133XVXTE9Dqefl33t4N56rV6/WjBkzlJOTo9zcXGVlZWn79u1av3696urqNGzYMF199dXt+dGc8fJn9Jvf/KZWrFihL774QnfeeacGDhyo8vJyrV69WtZaXXnllWEHn7GoLf7O1089i6cFAkJ5NaZjxozRSSedpM8++0x/+tOfNGDAAOXk5Ki4uDh4lubKK6/UUUcd1W6f7VAoajySnJysX//613r55Zc1d+5cLViwQOnp6Ro/fryuuOKKFl3klpWVpfvuu0/PP/+8FixYoPnz56tjx44655xzNGnSpLi4I67k7Zju2LFDa9asCeurqakJ69u7d69n2SORV+NZP3fe7/cHb2rWlHgoarwa0xNPPFFVVVVasWKF8vPztXv3biUnJ6t79+4666yzdO6558bsfSoO5OXfe3g3nqNGjdKOHTuUn5+vgoIC7dmzRx06dFBeXp7GjRunM844o9GFxLHKy5/R+n299tprmjNnjhYvXqykpCQNHz5c559/flxMsfb673xpaam++uorJSQk6JRTTmmj1JHNqzE1xuiuu+7S7Nmz9eGHH2rjxo3asGGD0tLSdNxxx+ncc8+NqBkExnLDDwAAAABRLD6+VgEAAAAQsyhqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENUoagAAAABENYoaAAAAAFGNogYAAABAVKOoAQAAABDVKGoAAAAARDWKGgAAAABRjaIGAAAAQFSjqAEAAAAQ1ShqAAAAAEQ1ihoAAAAAUY2iBgAAAEBUo6gBAAAAENX+P5FPjiuDPAgRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Value at both ends\n",
    "fa, fb = f(a), f(b)\n",
    "\n",
    "# Plot\n",
    "x = np.linspace(a, b, 81)\n",
    "\n",
    "plt.plot(x, f(x))\n",
    "plt.plot([a, b], [fa, fb])\n",
    "\n",
    "# Trapezoidal rule\n",
    "Itrap = 0.5*(fa + fb)*(b-a)\n",
    "print(Itrap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a69e5460",
   "metadata": {},
   "source": [
    "#### 다중 적분 (Composite trapezoid rule)\n",
    "- $[a, b]$ 구간을 여러 간격으로 (sub intervals) 나눠서 적분\n",
    "\n",
    ":::{figure-md} CompTrapezoid\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/d/d1/Integration_num_trapezes_notation.svg\" alt=\"Comptrapezoid-fig\">\n",
    "\n",
    "Composite Trapezoid rule (from Wikipedia)\n",
    ":::   \n",
    "\n",
    "- *(n+1)* 개 점으로 등간격 분할\n",
    "\n",
    "    * $h = (b-a)/n$\n",
    "    \n",
    "    $$\n",
    "    \\begin{align}\n",
    "    I &= \\int_{x_0}^{x_1} f(x) dx + \\int_{x_1}^{x_2} f(x) dx + ... + \\int_{x_{n-1}}^{x_n} f(x) dx \\\\\n",
    "    & \\approx h \\frac{f(x_0) + f(x_1)}{2} + h \\frac{f(x_1) + f(x_2)}{2} + ... \n",
    "    + h \\frac{f(x_{n-1}) + f(x_n)}{2} \\\\\n",
    "    &= \\frac{h}{2} \\left [\n",
    "    f(x_0) + 2\\sum_{i=1}^{n-1} f(x_i) + f(x_n)\n",
    "    \\right ]\n",
    "    \\end{align}\n",
    "    $$\n",
    "    \n",
    "    * 오차\n",
    "    \n",
    "    $$\n",
    "    E_t = \\sum_{i=1}^n -\\frac{1}{12} f''(\\xi_i) (h)^3 = -\\frac{(b-a)^3}{12n^3} \\sum_{i=1}^n f''(\\xi_i)\n",
    "    \\approx -\\frac{(b-a)^3}{12n^2} \\bar{f}''\n",
    "    $$\n",
    "    \n",
    "- 예제    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "105cfbb2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0688000000000115\n",
      "2, 0.4000, 1.0688, 34.9\n",
      "3, 0.2667, 1.3696, 16.5\n",
      "4, 0.2000, 1.4848, 9.5\n",
      "5, 0.1600, 1.5399, 6.1\n",
      "6, 0.1333, 1.5703, 4.3\n",
      "7, 0.1143, 1.5887, 3.2\n",
      "8, 0.1000, 1.6008, 2.4\n",
      "9, 0.0889, 1.6091, 1.9\n",
      "10, 0.0800, 1.6150, 1.6\n"
     ]
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "Ie = 1.640533\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Num of sub-interval\n",
    "n = 2\n",
    "\n",
    "# Trapezoid rule\n",
    "h = (b-a)/n\n",
    "I = (f(a) + 2*sum([f(a + i*h) for i in range(1, n)]) + f(b))* h/2\n",
    "print(I)\n",
    "\n",
    "# Increasing n\n",
    "for n in range(2, 11):\n",
    "    h = (b-a)/n\n",
    "    I = (f(a) + 2*sum([f(a + i*h) for i in range(1, n)]) + f(b))* h/2\n",
    "    et = abs((I-Ie)/Ie)*100\n",
    "    print(\"{}, {:.4f}, {:.4f}, {:.1f}\".format(n, h, I, et))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "804c1f69",
   "metadata": {},
   "source": [
    "#### DIY\n",
    "(Single and Composite) Trapezoid 적분 계산 함수를 구성하시오."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b111a17",
   "metadata": {},
   "source": [
    "### Simpson's Rule\n",
    "- 구간 내를 다항식으로 근사해서 적분 정확도 향상\n",
    "\n",
    "#### Simpson 1/3 rule\n",
    "- 구간 내에 3개의 등간격 점을 이용해서 2차 Lagrange 다항식으로 근사해서 적분\n",
    "\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    I &\\approx \\int_a^b f_2(x)dx \\\\\n",
    "    &= \\int_a^b \\left [\n",
    "    \\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} f(x_0)\n",
    "    +\\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} f(x_1)\n",
    "    +\\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} f(x_2)\n",
    "    \\right ] dx\\\\\n",
    "    &= \\frac{h}{3}[f(x_0) + 4f(x_1) + f(x_2)]\n",
    "    = \\frac{b-a}{6}[f(x_0) + 4f(x_1) + f(x_2)]\n",
    "    \\end{align}\n",
    "    $$\n",
    "    \n",
    "    * $h=(b-a)/2$\n",
    "\n",
    "\n",
    ":::{figure-md} Simpson\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/c/ca/Simpsons_method_illustration.svg\" alt=\"simpson-fig\">\n",
    "\n",
    "Simpson 1/3 rule (from Wikipedia)\n",
    ":::    \n",
    "\n",
    "- 오차\n",
    "    \n",
    "    $$\n",
    "    E_t = -\\frac{1}{90} h^5 f^{(4)}(\\xi)  = - \\frac{(b-a)^5}{2880} f^{(4)} (\\xi)\n",
    "    $$    \n",
    "    \n",
    "    * 유도\n",
    "        - Newton 보간식\n",
    "            - $\\alpha = \\frac{x-a}{h}$, $h=(b-a)/2$, $\\Delta^n f(a)$ 은 제차분식\n",
    "            \n",
    "            $$\\begin{align}\n",
    "            I &= \\int_0^2 \\left [\n",
    "            f(a) + \\Delta f(x_0)\\alpha + \\frac{\\Delta^2 f(x_0)}{2} \\alpha (\\alpha-1)\n",
    "            + \\frac{\\Delta^3 f(x_0)}{6} \\alpha (\\alpha-1)(\\alpha-2)\n",
    "            + \\frac{f^{(4)}(\\xi)}{24} \\alpha (\\alpha-1)(\\alpha-2)(\\alpha-3) h^4\n",
    "            \\right ] d \\alpha\n",
    "            \\end{align}\n",
    "            $$\n",
    "        \n",
    "            - $f^{(4)}(xi)$ 가 거의 변하지 않는 상수로 가정하면\n",
    "        \n",
    "            $$\\begin{align}                \n",
    "            I & \\approx h \\left [\n",
    "            \\alpha f(x_0) + \\frac{\\alpha^2}{2} \\Delta f(a) \n",
    "            + \\left (\n",
    "            \\frac{\\alpha^3}{6} - \\frac{\\alpha^2}{4}\n",
    "            \\right) \\Delta^2 f(x_0)\n",
    "            + \\left (\n",
    "            \\frac{\\alpha^4}{24} - \\frac{\\alpha^3}{6} + \\frac{\\alpha^2}{6}\n",
    "            \\right) \\Delta^3 f(x_0)\n",
    "            + \\left (\n",
    "            \\frac{\\alpha^5}{120}\n",
    "            - \\frac{\\alpha^4}{16} + \\frac{11\\alpha^3}{72} - \\frac{\\alpha^2}{8}\n",
    "            \\right) f^{(4)}(\\xi) h^4\n",
    "            \\right ]_0^2 \\\\\n",
    "            & = h \\left [\n",
    "            2f(a) + 2\\Delta f(x_0) + \\frac{\\Delta^2 f(x_0)}{3} + (0) \n",
    "            \\Delta^3 f(x_0) \n",
    "            \\right ] - \\frac{1}{90} f^{(4)}(\\xi) h^5 \\\\\n",
    "            & = \\frac{h}{3} [f(x_0) + 4 f(x_1) + f(b)]  - \\frac{1}{90} f^{(4)}(\\xi) h^5\n",
    "            \\end{align}\n",
    "            $$\n",
    "\n",
    "#### 다중 적분 (Composite  Simpson 1/3 rule)  \n",
    "\n",
    "- *(n+1)* 개 점으로 등간격 분할 (짝수)\n",
    "\n",
    "    * $h = (b-a)/n$\n",
    "    \n",
    "    $$\n",
    "    \\begin{align}\n",
    "    I &= \\int_{x_0}^{x_2} f(x) dx + \\int_{x_2}^{x_4} f(x) dx + ... + \\int_{x_{n-2}}^{x_n} f(x) dx \\\\\n",
    "    & \\approx 2h \\frac{f(x_0) + 4 f(x_1) + f(x_2)}{6} \n",
    "    + 2h \\frac{f(x_2) + 4 f(x_3) + f(x_4)}{6} + ...\n",
    "    + 2h \\frac{f(x_{n-2}) + 4 f(x_{n-1}) + f(x_n)}{6}\n",
    "    \\\\\n",
    "    &=\\frac{b-a}{3n} \\left [\n",
    "    f (x_0) + 4 \\sum_{i=1,3,5}^{n-1} f(x_i) + 2\\sum_{i=2,4,6}^{n-2} f(x_i) + f(x_n)\n",
    "    \\right ]\n",
    "    \\end{align}\n",
    "    $$\n",
    "    \n",
    "- 오차\n",
    "    \n",
    "    $$\n",
    "    E_a = -\\frac{(b-a)^5}{180 n^4} \\bar{f}^{(4)}\n",
    "    $$    \n",
    "    \n",
    "#### Simpson 3/8 rule\n",
    "- 구간 내에 4개의 등간격 점을 이용해서 2차 Lagrange 다항식으로 근사해서 적분\n",
    "\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    I &\\approx \\int_a^b f_3(x)dx \\\\\n",
    "    &=  \\frac{3h}{8}[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]\n",
    "    = \\frac{b-a}{8}[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]\n",
    "    \\end{align}\n",
    "    $$\n",
    "    \n",
    "    * $h=(b-a)/3$   \n",
    "\n",
    "- 오차\n",
    "    \n",
    "    $$\n",
    "    E_t = -\\frac{3}{80} h^5 f^{(4)}(\\xi)  = - \\frac{(b-a)^5}{6480} f^{(4)} (\\xi)\n",
    "    $$        \n",
    "    \n",
    "- 다중 적분시, 홀수개 분활에 사용 \n",
    "    \n",
    "#### 예제\n",
    "- 앞선 문제를 Simpson rule로 풀어보시오.\n",
    "\n",
    "- Composite Simpson 1/3 rule로 풀어보시오."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bc334a04",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3674666666666742\n",
      "1.519170370370378\n"
     ]
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Value at both ends\n",
    "fa, fb = f(a), f(b)\n",
    "\n",
    "# Simpson 1/3\n",
    "def simpson13(f, a, b):\n",
    "    x0, x1, x2 = a, 0.5*(a+b), b\n",
    "    \n",
    "    return (b-a)/6*(f(x0) + 4*f(x1) + f(x2))\n",
    "\n",
    "\n",
    "# Simpson 3/8\n",
    "def simpson38(f, a, b):\n",
    "    x0, x1, x2, x3 = a, (2*a + b)/3, (a +2*b)/3, b\n",
    "    \n",
    "    return (b-a)/8*(f(x0) + 3*f(x1) + 3*f(x2) + f(x3))\n",
    "\n",
    "\n",
    "I13 = simpson13(f, a, b)\n",
    "print(I13)\n",
    "\n",
    "I38 = simpson38(f, a, b)\n",
    "print(I38)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6f6b8b8c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6234666666666717"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "Ie = 1.640533\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Sub-interval\n",
    "n = 4\n",
    "xi = np.linspace(a, b, n + 1)\n",
    "\n",
    "# Composite integration\n",
    "(f(xi[0]) + 4*sum(f(xj) for xj in xi[1:-1:2]) + 2*sum(f(xj) for xj in xi[2:-2:2]) + f(xi[-1]))* (b-a)/ (3*n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "445afdbe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4, 1.6235, 1.0\n",
      "6, 1.6372, 0.2\n",
      "8, 1.6395, 0.1\n",
      "10, 1.6401, 0.0\n"
     ]
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "Ie = 1.640533\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "for n in [4, 6, 8, 10]:\n",
    "    xi = np.linspace(a, b, n + 1)\n",
    "    I = (f(x[0]) + 4*sum(f(xj) for xj in xi[1:-1:2]) + 2*sum(f(xj) for xj in xi[2:-2:2]) + f(xi[-1]))* (b-a)/ (3*n)\n",
    "    et = abs((I-Ie)/Ie)*100\n",
    "    print(\"{}, {:.4f}, {:.1f}\".format(n, I, et))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e011460",
   "metadata": {},
   "source": [
    "## Romberg Integration\n",
    "### Richardson Extrapolation\n",
    "- 정확도가 낮은 수치적분 기법을 활용하여 고차 수치 적분을 외삽\n",
    "\n",
    "- Trapezoid intgration을 사용하는 경우\n",
    "\n",
    "    $$\n",
    "        I_1 \\approx \\frac{h}{2} \\left [\n",
    "        f(x_0) + 2\\sum_{i=1}^{n-1} f(x_i) + f(x_n)\n",
    "        \\right ] + c_1 h^2 + c_2 h^4 + ...\n",
    "    $$\n",
    "\n",
    "    * 간격을 절반으로 하면\n",
    "    \n",
    "        $$\n",
    "        I_2 \\approx \\frac{h}{4} \\left [\n",
    "        f(x_0) + 2\\sum_{i=1}^{n-1} f(x_i) + f(x_n)\n",
    "        \\right ] + c_1 \\frac{h^2}{4} + c_2 \\frac{h^4}{16} + ...\n",
    "        $$\n",
    "        \n",
    "- 두 식에서 Leading error 항을 제거하기 위해\n",
    "\n",
    "    $$\n",
    "    I_{12} = \\frac{4 I_2 - I_1}{3} = I + \\frac{1}{4}c_2 h^4 + ...\n",
    "    $$\n",
    "    \n",
    "### Romberg integation\n",
    "- Rihcardson extrapolation 일반화\n",
    "\n",
    "    $$\n",
    "    I_{j,k} = \\frac{4^{k-1} I_{j+1, k-1} - I_{j,k-1}}{4^{k-l} - 1}\n",
    "    $$\n",
    "    \n",
    "#### 예제"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "772bbd76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2, 0.4000, 1.0688, 34.9\n",
      "4, 0.2000, 1.4848, 9.5\n",
      "8, 0.1000, 1.6008, 2.4\n",
      "[1.0688000000000115, 1.4848000000000068, 1.600800000000004]\n",
      "[1.623466666666672, 1.6394666666666697]\n",
      "1.6405333333333363\n"
     ]
    }
   ],
   "source": [
    "# Trapezoid rule with decreasing the size of sub-interval\n",
    "I1 = []\n",
    "for n in [2, 4, 8]:\n",
    "    h = (b-a)/n\n",
    "    I = (f(a) + 2*sum([f(a + i*h) for i in range(1, n)]) + f(b))* h/2\n",
    "    et = abs((I-Ie)/Ie)*100\n",
    "    print(\"{}, {:.4f}, {:.4f}, {:.1f}\".format(n, h, I, et))\n",
    "    I1.append(I)\n",
    "\n",
    "print(I1)        \n",
    "    \n",
    "# Romberg Integration\n",
    "k = 2\n",
    "I21 = (4**(k-1)*I1[1] -I1[0]) / (4**(k-1)-1)\n",
    "I22 = (4**(k-1)*I1[2] -I1[1]) / (4**(k-1)-1)\n",
    "I2 = [I21, I22]\n",
    "print(I2)\n",
    "\n",
    "# Romberg Integration\n",
    "k = 3\n",
    "I3 = (4**(k-1)*I2[1] -I2[0]) / (4**(k-1)-1)\n",
    "\n",
    "print(I3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad65d2b5",
   "metadata": {},
   "source": [
    "## Gauss (-Legendre) Integration\n",
    "- 적분 좌표와 Weighting 계수를 조절해여 고차 정확도 적분\n",
    "\n",
    "    $$\n",
    "    I = \\int_a^b f(x) dx = \\sum_i w_i f(\\xi_i)\n",
    "    $$\n",
    "    \n",
    "    * Transformation: $\\xi=\\frac{x - \\frac{a+b}{2}}{\\frac{b-a}{2}}$\n",
    "\n",
    "### 2 Points\n",
    "- 두점으로 3차 함수를 정확하게 적분하는 공식\n",
    "    * $f(x)=1, x, x^2, x^3$\n",
    "\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    w_0 f(x_0) + w_1 f(x_1) &= \\int_{-1}^1 1 dx = 2 \\\\\n",
    "    w_0 f(x_0) + w_1 f(x_1) &= \\int_{-1}^1 x dx = 0 \\\\\n",
    "    w_0 f(x_0) + w_1 f(x_1) &= \\int_{-1}^1 x^2 dx = \\frac{2}{3} \\\\\n",
    "    w_0 f(x_0) + w_1 f(x_1) &= \\int_{-1}^1 x^3 dx = 0 \n",
    "    \\end{align}\n",
    "    $$\n",
    "\n",
    "- 연립방정식의 해\n",
    "    * $w_0 = w_1 = 1$\n",
    "    * $x_0 = -\\frac{1}{\\sqrt{3}}, x_1 = \\frac{1}{\\sqrt{3}}$\n",
    "    \n",
    "- 일반 구간\n",
    "\n",
    "    $$\n",
    "    I = \\int_a^b f(x) dx = \\frac{b-a}{2} \\int_{-1}^{1} f(\\xi) d\\xi\n",
    "    $$\n",
    "\n",
    "### 다점 Gauss-Legendre Quadrature rule\n",
    "- See [Wikipeida](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature)\n",
    "\n",
    "#### 예제"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "98278e46",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8225777777777772\n"
     ]
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "Ie = 1.640533\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Transformation [-1, 1 ] -> [a, b]\n",
    "xi_to_x = lambda xi: (a+b)/2 + xi*(b-a)/2\n",
    "\n",
    "# 2 Points rule\n",
    "w0=w1 = 1\n",
    "xi0, xi1 = -1/np.sqrt(3), 1/np.sqrt(3)\n",
    "\n",
    "# Gauss Qaudrature\n",
    "I = (w0*f(xi_to_x(xi0)) + w1*f(xi_to_x(xi1))) *(b-a) /2\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1af50e5a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6405333333333294\n"
     ]
    }
   ],
   "source": [
    "# function\n",
    "f = lambda x: 0.2 + 25*x - 200*x**2 + 675*x**3 -900*x**4 +400*x**5\n",
    "Ie = 1.640533\n",
    "\n",
    "# Interval\n",
    "a, b = 0, 0.8\n",
    "\n",
    "# Transformation [-1, 1 ] -> [a, b]\n",
    "xi_to_x = lambda xi: (a+b)/2 + xi*(b-a)/2\n",
    "\n",
    "# 3 points rule\n",
    "xi = -np.sqrt(3/5), 0 , np.sqrt(3/5)\n",
    "w = 5/9, 8/9, 5/9\n",
    "\n",
    "# Gauss quadrature\n",
    "I = 0\n",
    "for wi, xii in zip(w, xi):\n",
    "    I += wi*f(xi_to_x(xii))\n",
    "    \n",
    "I *= (b - a)/2\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f86d818e",
   "metadata": {},
   "source": [
    "## Scipy 활용\n",
    "\n",
    "Scipy 내 `integrate` 모둘에서 수치 적분에 대한 여러 함수를 제공한다.\n",
    "\n",
    "- 수치 적분: `scipy.integrate` 모듈\n",
    "    * https://docs.scipy.org/doc/scipy/tutorial/integrate.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}